⚙️ How-To Coastal & marine Geometry & topology Coastal & Marine · step 12 of 20 Updated

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. A country such as Fiji is stored as a single feature made of 44 parts spanning the full −180 to 180 range, with 18 parts west of −170° and 26 east of 170°. Its naive centroid longitude is 163.41°E, which is over a thousand kilometres from the country, and every clip, buffer and bounding-box operation inherits that.

Quick answer

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

import geopandas as gpd, numpy as np

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))

clipped = gpd.clip(data.to_crs(zone.crs), zone)          # cuts geometries at the boundary
inside = gpd.sjoin(data.to_crs(zone.crs), zone, predicate="within")   # keeps whole features

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.

Three scenes showing clip cutting a feature, within dropping it and intersects keeping it whole.
The same dataset and the same zone give three different results.

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.

Two panels showing a Fiji-like territory split at the antimeridian and the same territory in a 0 to 360 frame.
The same geometry, two longitude frames; every operation works in the second.

Code examples

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(after / before, 4) if before else None}
    return out, report

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

zone = gpd.read_file("eez.gpkg").query("territory == 'Fiji'")
print("naive centroid longitude:", round(float(zone.geometry.centroid.x.iloc[0]), 2))

z360 = shift_longitudes(zone)
d360 = shift_longitudes(data.to_crs(zone.crs))
clipped = gpd.clip(d360, z360)
result = shift_longitudes(clipped, to_360=False)
print("centroid in the 0–360 frame:", round(float(z360.geometry.centroid.x.iloc[0]), 2))
naive centroid longitude: 163.41
centroid in the 0–360 frame: 178.52

163.41°E is in the open Pacific, roughly 1,500 km west of Fiji. 178.52° is the country, whose shifted bounds run from 174.59 to 181.78. The shifted frame is what every operation needs.

Example 3 — area statistics per zone, correctly projected 

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

def per_zone_area(features, zones, name_col="territory"):
    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
        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],
            "features": len(clipped),
            "zone_km2": round(g.to_crs(crs).area.iloc[0] / 1e6),
            "feature_km2": round(clipped.to_crs(crs).area.sum() / 1e6, 1) if len(clipped) else 0.0,
        })
    return gpd.pd.DataFrame(rows)

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.

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 their average is −1 — the Atlantic. Fiji's 44 parts average to 163.41, which is not as dramatic and is still 1,500 km from the country. 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.

Checklist of clip checks: matching CRS, antimeridian handling, predicate choice, equal-area projection, overlapping zones and attribute validity after cutting.
A cut polygon keeps the whole feature’s attributes, which is rarely what you want.

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.

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.