How to clip a dataset to an exclusive economic zone

Problem statement

Clipping to an EEZ is a gpd.clip call, and three things around it decide whether the result is right: which EEZ layer, whether the zone crosses the antimeridian, and whether you need the features that straddle the boundary kept, cut or dropped.

The antimeridian is the one that fails catastrophically. In the Marine Regions EEZ layer, Fiji's zone is a single feature of two parts cut at 180°: one from 172.76°E to 180°, the other from −180° to −176.27°. Its bounding box therefore spans the whole −180 to 180 range, and its naive centroid longitude is 93.68°E — in the Indian Ocean. The clip itself survives, because each part is a valid polygon on its own side of the line. The bounding box, the centroid and everything built on them do not.

Quick answer

Check for an antimeridian crossing before anything else, then clip with the right predicate:

import geopandas as gpd


def crosses_antimeridian(gdf, margin=10.0):
    b = gdf.explode(index_parts=False).bounds
    return bool((b.minx < -180 + margin).any() and (b.maxx > 180 - margin).any())


eez = gpd.read_file("eez.gpkg")
zone = eez[eez["territory"] == "Fiji"]
print("crosses the antimeridian:", crosses_antimeridian(zone))

data = gpd.read_file("data.gpkg")                        # M4.5+ earthquakes, 2015-2024
print("CRS:", data.crs.to_epsg(), "against", zone.crs.to_epsg())
data = data.to_crs(zone.crs)
clipped = gpd.clip(data, zone)                           # cuts geometries at the boundary
inside = gpd.sjoin(data, zone, predicate="within")       # keeps whole features
print(f"{len(data):,} earthquakes, {len(clipped):,} clipped, {len(inside):,} within")

crosses the antimeridian: True
CRS: 4979 against 4326
13,512 earthquakes, 2,538 clipped, 2,538 within

The earthquake catalogue comes as EPSG:4979, latitude–longitude with a height axis. geopandas warns about the mismatch rather than failing; convert it and the warning goes.

clip cuts; sjoin with within keeps only features entirely inside; sjoin with intersects keeps any feature that touches. Those are three different answers and the choice belongs to the analysis.

Step-by-step solution

1. Choose the zone layer deliberately

An authoritative EEZ layer from Marine Regions or a national hydrographic office carries the treaty boundaries and the disputed areas. A 200 NM buffer computed from a coastline is an approximation that understates the zone wherever straight baselines apply. Say which you used.

2. Detect an antimeridian crossing

A geometry whose parts reach both −180 and +180 crosses it. So does a bounding box that spans nearly the whole world for a country that does not.

3. Shift to a 0–360 frame if it does

The simplest reliable fix is to translate every negative longitude by +360, do the work, and translate back. Both the zone and the data have to be shifted, and the result has to be shifted back before writing.

from shapely.ops import transform

def to_0_360(geom):
    return transform(lambda x, y, z=None: (np.where(x < 0, x + 360, x), y), geom)

4. Decide what happens to straddling features

• Cut them — gpd.clip. Right for area-based analysis; it invents boundaries in the attribute table, so any per-feature attribute is no longer valid for the cut piece.
• Keep whole features that are entirely inside — predicate="within". Right when the features are indivisible.
• Keep whole features that touch — predicate="intersects". Right for a screening filter.

5. Reproject before area work, not after

Areas and distances computed in degrees are meaningless. Clip in a geographic CRS if you like, then reproject to an equal-area projection for any statistic.

6. Watch the overlaps

EEZs overlap where entitlements conflict, and a point can legitimately be in two. A spatial join then duplicates rows, which is correct and surprising.

7. Check the counts and the area against the zone

Features dropped, features cut, area retained. A clip that keeps 3% of the data usually means a CRS mismatch rather than a small zone.

Code examples

The data used on this page

Seven South Pacific zones from the Marine Regions EEZ layer and its Longhurst biogeographic provinces, both through the Marine Regions WFS, and every M4.5+ earthquake from the USGS catalogue between 5° and 30° S and 160°E and 160°W for 2015–2024:

import io

import geopandas as gpd
import pandas as pd
import requests

WFS = "https://geo.vliz.be/geoserver/MarineRegions/wfs"


def wfs(layer, where=None):
    params = {"service": "WFS", "version": "1.0.0", "request": "GetFeature",
              "typeName": f"MarineRegions:{layer}", "outputFormat": "application/json"}
    if where:
        params["CQL_FILTER"] = where
    return gpd.read_file(io.BytesIO(requests.get(WFS, params=params, timeout=300).content))


names = ["Fiji", "Tonga", "Vanuatu", "Samoa", "Wallis and Futuna", "New Caledonia", "Tuvalu"]
eez = pd.concat(wfs("eez", f"territory1='{n}'") for n in names)
eez = eez[["mrgid", "geoname", "territory1", "area_km2", "geometry"]].rename(columns={"territory1": "territory"})
eez.to_file("eez.gpkg")
wfs("longhurst")[["provcode", "provdescr", "geometry"]].to_file("longhurst.gpkg")

quakes = requests.get("https://earthquake.usgs.gov/fdsnws/event/1/query", timeout=300, params={
    "format": "geojson", "starttime": "2015-01-01", "endtime": "2025-01-01", "minmagnitude": 4.5,
    "minlatitude": -30, "maxlatitude": -5, "minlongitude": 160, "maxlongitude": 200})
gpd.read_file(io.BytesIO(quakes.content))[["id", "mag", "time", "geometry"]].to_file("data.gpkg")
print(eez[["territory", "area_km2"]].to_string(index=False))

        territory  area_km2
             Fiji   1302114
            Tonga    666052
          Vanuatu    625850
            Samoa    130480
Wallis and Futuna    250614
    New Caledonia   1175971
           Tuvalu    753133

The USGS query takes longitudes past 180 and returns them wrapped to −180–180, so the catalogue itself arrives split across the line.

Example 1 — a clip that reports what it did

import geopandas as gpd, numpy as np

def clip_to_zone(data, zone, mode="clip", area_crs=None):
    if data.crs != zone.crs:
        data = data.to_crs(zone.crs)
    z = zone.union_all()

    if mode == "clip":
        out = gpd.clip(data, zone)
        cut = int((~out.geometry.geom_equals(
            data.loc[out.index, "geometry"], align=True)).sum())
    elif mode == "within":
        out = data[data.within(z)]
        cut = 0
    else:
        out = data[data.intersects(z)]
        cut = 0

    report = {"mode": mode, "in": len(data), "out": len(out),
              "dropped": len(data) - len(out), "cut": cut}
    if area_crs is not None and (data.geom_type.isin(["Polygon", "MultiPolygon"])).all():
        before = data.to_crs(area_crs).area.sum() / 1e6
        after = out.to_crs(area_crs).area.sum() / 1e6
        report |= {"area_before_km2": round(before), "area_after_km2": round(after),
                   "area_retained": round(float(after / before), 4) if before else None}
    return out, report


provinces = gpd.read_file("longhurst.gpkg")
fiji_laea = "+proj=laea +lat_0=-18 +lon_0=178 +datum=WGS84 +units=m"
for mode in ("clip", "within", "intersects"):
    out, report = clip_to_zone(provinces, zone, mode=mode, area_crs=fiji_laea)
    print(report, sorted(out.provcode))

{'mode': 'clip', 'in': 54, 'out': 2, 'dropped': 52, 'cut': 2, 'area_before_km2': 362034913, 'area_after_km2': 1293903, 'area_retained': 0.0036} ['SPSG', 'WARM']
{'mode': 'within', 'in': 54, 'out': 0, 'dropped': 54, 'cut': 0, 'area_before_km2': 362034913, 'area_after_km2': 0, 'area_retained': 0.0} []
{'mode': 'intersects', 'in': 54, 'out': 2, 'dropped': 52, 'cut': 0, 'area_before_km2': 362034913, 'area_after_km2': 48472481, 'area_retained': 0.1339} ['SPSG', 'WARM']

Two Longhurst provinces reach Fiji's zone — SPSG (South Pacific Subtropical Gyre) and WARM (Western Pacific Warm Pool) — and neither lies wholly inside it, so within returns nothing. clip keeps 1,293,903 km², close to the zone's own area; intersects keeps both provinces whole, 48.5 million km², 37 times the zone. The cut pieces still carry each province's full description, which is harmless here and wrong for any attribute that scales with area.

Example 2 — the antimeridian, handled

import numpy as np, geopandas as gpd
from shapely.ops import transform

def shift_longitudes(gdf, to_360=True):
    def fn(x, y, z=None):
        x = np.asarray(x)
        return (np.where(x < 0, x + 360, x) if to_360
                else np.where(x > 180, x - 360, x), y)
    out = gdf.copy()
    out["geometry"] = out.geometry.apply(lambda g: transform(fn, g))
    return out

from pyproj import Geod

zone = gpd.read_file("eez.gpkg").query("territory == 'Fiji'")
naive = zone.geometry.iloc[0].centroid                   # geographic CRS: geopandas warns
z360 = shift_longitudes(zone)
d360 = shift_longitudes(data)
shifted = z360.geometry.iloc[0].centroid
off_km = Geod(ellps="WGS84").inv(naive.x, naive.y, shifted.x, shifted.y)[2] / 1000
print(f"naive centroid        {naive.x:7.2f}, {naive.y:6.2f}")
print(f"0-360 centroid        {shifted.x:7.2f}, {shifted.y:6.2f}   ({off_km:,.0f} km apart)")
print("bounds, -180 frame   ", zone.total_bounds.round(2))
print("bounds, 0-360 frame  ", z360.total_bounds.round(2))

for label, z, d in [("-180 frame", zone, data), ("0-360 frame", z360, d360)]:
    x0, y0, x1, y1 = z.total_bounds
    in_box = d.cx[x0:x1, y0:y1]
    print(f"{label:12} in bounding box {len(in_box):6,}   clipped {len(gpd.clip(d, z)):,}")

result = shift_longitudes(gpd.clip(d360, z360), to_360=False)
print("back in -180-180:", result.total_bounds.round(2))

naive centroid          93.68, -17.95
0-360 centroid         177.69, -17.95   (8,806 km apart)
bounds, -180 frame    [-180.    -25.1   180.     -9.78]
bounds, 0-360 frame   [172.76 -25.1  183.73  -9.78]
-180 frame   in bounding box 11,512   clipped 2,538
0-360 frame  in bounding box  3,311   clipped 2,538
back in -180-180: [-180.    -24.35  180.    -13.04]

The naive centroid is in the Indian Ocean. In the 0–360 frame it lands on Viti Levu, Fiji's main island — a hole in the zone, since EEZ polygons leave land out — and the bounding box shrinks from the whole planet to the zone, so a box pre-filter that let through 85% of a regional catalogue lets through only the earthquakes near Fiji. The clip count is the same in both frames — the thing that broke was never gpd.clip.

Example 3 — area statistics per zone, correctly projected

import geopandas as gpd, numpy as np
from pyproj import CRS

def zone_centre(g):
    """Centroid of a one-row GeoDataFrame, taken in the 0-360 frame when it crosses 180°."""
    if crosses_antimeridian(g):
        c = shift_longitudes(g).geometry.iloc[0].centroid
        return (c.x + 180) % 360 - 180, c.y
    c = g.geometry.iloc[0].centroid
    return c.x, c.y


def per_zone_area(features, zones, name_col="territory", naive=False):
    rows = []
    for _, zone in zones.iterrows():
        g = gpd.GeoDataFrame(geometry=[zone.geometry], crs=zones.crs)
        cx, cy = (zone.geometry.centroid.x, zone.geometry.centroid.y) if naive else zone_centre(g)
        crs = CRS.from_proj4(f"+proj=laea +lat_0={cy:.4f} +lon_0={cx:.4f} "
                             "+datum=WGS84 +units=m")
        clipped = gpd.clip(features.to_crs(zones.crs), g)
        rows.append({
            name_col: zone[name_col],
            "centre_lon": round(cx, 2),
            "features": len(clipped),
            "zone_km2": round(g.to_crs(crs).area.iloc[0] / 1e6),
            "listed_km2": zone["area_km2"],
            "feature_km2": round(clipped.to_crs(crs).area.sum() / 1e6) if len(clipped) else 0,
        })
    return gpd.pd.DataFrame(rows)


print(per_zone_area(provinces, eez).to_string(index=False))
print(per_zone_area(provinces, eez.query("territory == 'Fiji'"), naive=True).to_string(index=False))

        territory  centre_lon  features  zone_km2  listed_km2  feature_km2
             Fiji      177.69         2   1301938     1302114      1293912
            Tonga     -174.77         2    666049      666052       665928
          Vanuatu      168.56         3    625849      625850       624347
            Samoa     -172.75         1    130479      130480       130255
Wallis and Futuna     -176.81         2    250613      250614       250521
    New Caledonia      163.20         2   1175917     1175971      1175172
           Tuvalu      178.18         1    753131      753133       751841
territory  centre_lon  features  zone_km2  listed_km2  feature_km2
     Fiji       93.68         2   1305512     1302114      1297546

A per-zone Lambert azimuthal equal-area projection centred on each zone is the honest way to compare areas across a set of EEZs — a single global equal-area projection distorts shape badly at high latitudes even though the areas are right. The centre has to come from the 0–360 frame, though. Centred on the naive centroid at 93.68°E, Fiji's zone measured 1,305,512 km², 3,398 km² more than Marine Regions lists; centred at 177.69° it measured 1,301,938 km², within 176 km². An equal-area projection keeps areas exact, but the zone's edges are straight lines between vertices, and 84° from the centre those lines enclose a different shape.

Explanation

Why a naive centroid lands in the wrong ocean

The centroid is computed from the coordinates as numbers. A territory with parts at 177°E and 179°W is stored as parts at +177 and −179, and the area-weighted average falls somewhere between them — anywhere from the Atlantic to the western Pacific, depending on how the area splits. Three quarters of Fiji's zone lies between 172.76°E and 180° and a quarter between −180° and −176.27°; weighted by area, 176 and −178 average to 93.68°E, 8,806 km from the 0–360 centroid. Nothing about the value is flagged as impossible.

Why shifting beats splitting

You can split every geometry at the antimeridian and work with the pieces, and some pipelines do. Shifting to 0–360 is simpler, keeps each feature whole, and works for every operation — buffers, centroids, clips, bounding boxes — provided everything in the computation is shifted the same way and shifted back before writing. The only trap is forgetting one of the layers.

Why clip and sjoin give different counts

clip intersects the geometries, so a polygon straddling the boundary appears with a smaller area and its original attributes — which is fine for area and wrong for anything per-feature, such as a population count. within drops it. intersects keeps it whole, including the part outside the zone. Each is right for a different question and none is a default.

Why overlaps duplicate rows

Where two EEZs overlap, a point inside both matches both in a spatial join, and the join emits a row per match. That is the correct result and it breaks any downstream code that assumes one row per input feature. Either resolve the assignment explicitly or aggregate with a rule.

Edge cases or notes

• Check the CRS on both layers. A near-empty clip is usually a mismatch.
• clip keeps the attributes of a cut feature. They may no longer be true.
• Zones can be multipart with islands thousands of kilometres apart.
• Some EEZs are disputed and appear in more than one layer's polygons.
• Buffer-derived zones understate where straight baselines apply.
• Reproject for area, always. Degrees-squared is not a unit.
• Write in −180 to 180 unless the consumer expects otherwise.
• Record the zone source and version with the output.

Internal links

• A coastline or EEZ breaks at the antimeridian — the failure in detail
• Maritime boundaries explained: EEZ, territorial sea and baselines — where the zone comes from
• How to clip spatial data in Python with GeoPandas — the general operation
• How to select features by location in GeoPandas — the predicate choice
• Longitude conventions explained — 0–360 versus −180 to 180
• How to choose a projected CRS for an area — picking the equal-area projection
• Spatial join cardinality explained — why overlaps duplicate rows
• How to calculate area and distance in GeoPandas — measuring the result

FAQ

How do I clip data to an EEZ?

gpd.clip(data, zone) after matching the CRS — but check first whether the zone crosses the antimeridian, and decide whether straddling features should be cut, dropped or kept whole.

Why is my Pacific EEZ clip empty or enormous?

Because the zone crosses the antimeridian. Shift every layer to a 0–360 longitude frame, do the work, and shift back.

Should I use clip or a spatial join?

clip cuts geometries at the boundary, within keeps only whole features inside, intersects keeps whole features that touch. Pick by whether your features are divisible.

Why does my spatial join return more rows than it started with?

Because EEZs can overlap where entitlements conflict, and a feature inside two zones matches both.

Can I just buffer the coastline by 200 NM?

For an approximation, yes, in a suitable projection. It understates the zone wherever straight baselines apply and is not an authoritative boundary.

How do I compute areas per zone?

Reproject each zone and its clipped features to a Lambert azimuthal equal-area projection centred on that zone.