Skip to content
GeoJason

Spatial Operations Explained

Learn buffer clip union intersection and difference operations with real GeoJSON examples you can paste into tools.

What Are Spatial Operations?

Spatial operations transform, combine, or analyze geometries based on their geographic relationships. Attribute queries filter by property values; spatial operations work directly with coordinate geometry, computing intersections, measuring distances, merging shapes, or deriving new geometries from existing ones. The six operations below cover the majority of real-world GeoJSON processing workflows.

Each one is shown here in four languages: Python with Shapely, PostGIS, GDAL on the command line, and Turf.js in Node. Every example reads a.geojson, and b.geojson where the operation takes two inputs, and writes GeoJSON to standard output, so you can copy one and run it against your own files.

The output under each example is what that code printed. All 24 samples are executed against the inputs shown, and the page is built from the results. Where the four agree, they agree exactly. Where they disagree, the output under each tab is different, and that is worth knowing before you assume the four are interchangeable.

Buffer

A buffer creates a polygon at a fixed distance around a geometry. Input any point, line, or polygon; specify a radius; receive a new polygon covering the area within that distance. Buffers are the foundation of proximity analysis.

Real Scenario: 1km Exclusion Zone Around a School

A city planning department needs to identify all liquor license applications within 1km of any school. The input is a Point feature at the school's location. The output is a polygon approximating a circle with a radius of 1,000 meters. Any application whose point falls inside that polygon is flagged for review.

Input: School location (a.geojson):

json{
  "type": "Feature",
  "properties": {
    "name": "PS 116 Mary Lindley Murray",
    "type": "elementary_school"
  },
  "geometry": {
    "type": "Point",
    "coordinates": [
      -73.9857,
      40.7484
    ]
  }
}
python · Shapely 2.ximport json
from shapely.geometry import shape, mapping
from shapely.ops import transform
from pyproj import Transformer

school = json.load(open("a.geojson"))
point = shape(school["geometry"])

# Shapely works in the units it is given, and these units are degrees. A
# 1000 "degree" buffer is meaningless, so project to a metre-based CRS
# centred on the point, buffer there, and come back.
lon, lat = point.x, point.y
local = f"+proj=aeqd +lat_0={lat} +lon_0={lon} +datum=WGS84 +units=m"
to_m = Transformer.from_crs("EPSG:4326", local, always_xy=True).transform
to_deg = Transformer.from_crs(local, "EPSG:4326", always_xy=True).transform

ring = transform(to_deg, transform(to_m, point).buffer(1000, quad_segs=16))

print(json.dumps({
    "type": "Feature",
    "properties": {"source": school["properties"]["name"], "radius_m": 1000},
    "geometry": mapping(ring),
}))

What it printed

{
  "type": "Feature",
  "properties": {
    "source": "PS 116 Mary Lindley Murray",
    "radius_m": 1000
  },
  "geometry": {
    "type": "Polygon",
    "coordinates": [
      [
        [
          -73.97385927762036,
          40.74839939263463
        ],
        [
          -73.97391644973541,
          40.74751675113071
        ],
        [
          -73.97408709938702,
          40.74664262131516
        ],
        [
          -73.97436957881216,
          40.74578542104346
        ],
        [
          -73.97476116360465,
          40.74495340482346
        ],
        [
          -73.97525807910547,
          40.74415458435088
        ],
        [
          -73.97585553688853,
          40.74339665139346
        ],
        [
          -73.97654778098554,
          40.74268690376545
        ],
        [
          -73.97732814340061,
          40.742032175103134
        ],
        [
          -73.97818910837695,
          40.741438769115426
        ],
        [
          -73.97912238479437,
          40.7409123989405
        ],
        [
          -73.98011898600032,
          40.740458132189886
        ],
        [
          -73.9811693163059,
          40.74008034220692
        ],
        [
          -73.98226326331584,
          40.7397826660064
        ],
        [
          -73.98339029520659,
          40.739567969298385
        ],
        [
          -73.98453956201968,
          40.73943831893117
        ],
        [
          -73.9857,
          40.73939496301705
        ],
        [
          -73.98686043798033,
          40.73943831893117
        ],
        [
          -73.9880097047934,
          40.739567969298385
        ],
        [
          -73.98913673668416,
          40.7397826660064
        ],
        [
          -73.99023068369411,
          40.74008034220692
        ],
        [
          -73.99128101399967,
          40.740458132189886
        ],
        [
          -73.99227761520564,
          40.740912398940495
        ],
        [
          -73.99321089162305,
          40.741438769115426
        ],
        [
          -73.99407185659939,
          40.742032175103134
        ],
        [
          -73.99485221901446,
          40.74268690376545
        ],
        [
          -73.99554446311146,
          40.74339665139346
        ],
        [
          -73.9961419208945,
          40.74415458435089
        ],
        [
          -73.99663883639533,
          40.74495340482346
        ],
        [
          -73.99703042118783,
          40.74578542104346
        ],
        [
          -73.99731290061297,
          40.74664262131516
        ],
        [
          -73.99748355026459,
          40.74751675113072
        ],
        [
          -73.99754072237964,
          40.74839939263463
        ],
        [
          -73.99748386187272,
          40.749282045673475
        ],
        [
          -73.9973135118543,
          40.75015620965054
        ],
        [
          -73.99703130857272,
          40.75101346539755
        ],
        [
          -73.99663996582207,
          40.7518455562747
        ],
        [
          -73.99614324895984,
          40.752644467717325
        ],
        [
          -73.99554593877859,
          40.7534025044617
        ],
        [
          -73.99485378557435,
          40.75411236470513
        ],
        [
          -73.99407345385,
          40.75476721048357
        ],
        [
          -73.99321245818295,
          40.75536073358742
        ],
        [
          -73.9922790908728,
          40.75588721637779
        ],
        [
          -73.99128234206505,
          40.756341586915354
        ],
        [
          -73.99023181312087,
          40.75671946786835
        ],
        [
          -73.98913762406907,
          40.75701721872607
        ],
        [
          -73.98801031603477,
          40.757231970909366
        ],
        [
          -73.98686074958847,
          40.757361655438096
        ],
        [
          -73.9857,
          40.75740502288715
        ],
        [
          -73.98453925041152,
          40.757361655438096
        ],
        [
          -73.98338968396523,
          40.757231970909366
        ],
        [
          -73.98226237593092,
          40.75701721872607
        ],
        [
          -73.98116818687913,
          40.75671946786835
        ],
        [
          -73.98011765793494,
          40.75634158691535
        ],
        [
          -73.97912090912719,
          40.75588721637778
        ],
        [
          -73.97818754181704,
          40.75536073358742
        ],
        [
          -73.97732654615,
          40.75476721048357
        ],
        [
          -73.97654621442564,
          40.75411236470513
        ],
        [
          -73.9758540612214,
          40.753402504461704
        ],
        [
          -73.97525675104015,
          40.752644467717325
        ],
        [
          -73.97476003417793,
          40.7518455562747
        ],
        [
          -73.97436869142729,
          40.75101346539754
        ],
        [
          -73.9740864881457,
          40.75015620965054
        ],
        [
          -73.97391613812727,
          40.749282045673475
        ],
        [
          -73.97385927762036,
          40.74839939263463
        ]
      ]
    ]
  }
}

1 feature, Polygon, 65 positions, 3.137 km²

This is the one operation where the four disagree. Each returns a 65-position polygon and each covers a different area, from 3,136,548 m² to 3,140,416 m², a spread of 0.12% on what is nominally the same 1km circle.

The reason is that a circle on a sphere is not a circle on a plane, and each library approximates it differently. Shapely and GDAL are both projecting to a local azimuthal equidistant coordinate system and buffering there, so they agree to the square metre. PostGIS computes on the spheroid through its geography type. Turf works on the sphere directly. None is wrong; they answer slightly different questions, and the difference is about four thousand square metres on a three square kilometre circle.

Two practical consequences. A buffer distance in degrees is meaningless, because a degree of longitude is about 111 km at the equator and near zero at the poles, which is why three of the four examples above have a projection step in them. And if you buffer the same data in two tools and compare the results, expect them to differ. Use the GeoJSON buffer tool to generate buffers with a configurable radius.

Clip

A clip, also called a cookie-cutter or mask operation, trims a set of input features to the boundary of a clipping polygon. Features entirely inside the boundary pass through unchanged. Features entirely outside are removed. Features that cross the boundary are split, keeping only the portion inside.

Real Scenario: Extract Features Within a City Boundary

A transit agency wants to extract bike lanes from a regional dataset to produce a city-specific layer. The input is a LineString that runs out of the city on both sides, plus a Polygon for the city boundary. The output is the middle section, cut exactly at the boundary.

Input 1: Bike lane crossing the city boundary (a.geojson):

json{
  "type": "Feature",
  "properties": {
    "name": "Coastal Trail Connector",
    "surface": "asphalt"
  },
  "geometry": {
    "type": "LineString",
    "coordinates": [
      [
        -122.49,
        37.8199
      ],
      [
        -122.4783,
        37.8199
      ],
      [
        -122.465,
        37.8199
      ]
    ]
  }
}

Input 2: City boundary (b.geojson):

json{
  "type": "Feature",
  "properties": {
    "name": "City Boundary"
  },
  "geometry": {
    "type": "Polygon",
    "coordinates": [
      [
        [
          -122.485,
          37.81
        ],
        [
          -122.47,
          37.81
        ],
        [
          -122.47,
          37.828
        ],
        [
          -122.485,
          37.828
        ],
        [
          -122.485,
          37.81
        ]
      ]
    ]
  }
}
python · Shapely 2.ximport json
from shapely.geometry import shape, mapping

line = json.load(open("a.geojson"))
boundary = json.load(open("b.geojson"))

clipped = shape(line["geometry"]).intersection(shape(boundary["geometry"]))

print(json.dumps({
    "type": "Feature",
    "properties": line["properties"],
    "geometry": mapping(clipped),
}))

What it printed

{
  "type": "Feature",
  "properties": {
    "name": "Coastal Trail Connector",
    "surface": "asphalt"
  },
  "geometry": {
    "type": "LineString",
    "coordinates": [
      [
        -122.485,
        37.8199
      ],
      [
        -122.4783,
        37.8199
      ],
      [
        -122.47,
        37.8199
      ]
    ]
  }
}

1 feature, LineString, 3 positions, 1,321 m long

All four return the same line, cut at the boundary at both ends, 1,321 metres long. The two end positions are new; the middle one was already there.

One caveat on the JavaScript example. Turf has no line-on-polygon clip, so it clips to the boundary's bounding box instead. That gives the same answer here because this boundary is a rectangle. Against any other shape it keeps parts of the line that fall outside the polygon but inside its box, which is the wrong answer and a silent one. Reach for Shapely, PostGIS or GDAL when the boundary is irregular.

A line that leaves the boundary and comes back returns as a MultiLineString, so handle both types. The clip keeps the original properties. Use the GeoJSON clip tool to clip a FeatureCollection to any polygon boundary.

Union

Union merges two or more overlapping or adjacent polygons into a single polygon covering the combined area. Boundaries between the inputs dissolve where they touch or overlap, and the outer perimeter of the combined shape becomes the output geometry.

Real Scenario: Merge Overlapping Land Parcels

A developer has acquired two adjacent parcels that partially overlap because of a surveying discrepancy. To calculate the total buildable area and produce one parcel boundary for permit applications, they union the two polygons.

Input 1: Parcel A-1042 (a.geojson):

json{
  "type": "Feature",
  "properties": {
    "parcel_id": "A-1042"
  },
  "geometry": {
    "type": "Polygon",
    "coordinates": [
      [
        [
          -73.987,
          40.748
        ],
        [
          -73.985,
          40.748
        ],
        [
          -73.985,
          40.7495
        ],
        [
          -73.987,
          40.7495
        ],
        [
          -73.987,
          40.748
        ]
      ]
    ]
  }
}

Input 2: Parcel A-1043 (b.geojson):

json{
  "type": "Feature",
  "properties": {
    "parcel_id": "A-1043"
  },
  "geometry": {
    "type": "Polygon",
    "coordinates": [
      [
        [
          -73.9855,
          40.7488
        ],
        [
          -73.9835,
          40.7488
        ],
        [
          -73.9835,
          40.7503
        ],
        [
          -73.9855,
          40.7503
        ],
        [
          -73.9855,
          40.7488
        ]
      ]
    ]
  }
}
python · Shapely 2.ximport json
from shapely.geometry import shape, mapping
from shapely.ops import unary_union

a = json.load(open("a.geojson"))
b = json.load(open("b.geojson"))

merged = unary_union([shape(a["geometry"]), shape(b["geometry"])])

print(json.dumps({
    "type": "Feature",
    "properties": {"parcel_ids": f"{a['properties']['parcel_id']},{b['properties']['parcel_id']}"},
    "geometry": mapping(merged),
}))

What it printed

{
  "type": "Feature",
  "properties": {
    "parcel_ids": "A-1042,A-1043"
  },
  "geometry": {
    "type": "Polygon",
    "coordinates": [
      [
        [
          -73.985,
          40.748
        ],
        [
          -73.987,
          40.748
        ],
        [
          -73.987,
          40.7495
        ],
        [
          -73.9855,
          40.7495
        ],
        [
          -73.9855,
          40.7503
        ],
        [
          -73.9835,
          40.7503
        ],
        [
          -73.9835,
          40.7488
        ],
        [
          -73.985,
          40.7488
        ],
        [
          -73.985,
          40.748
        ]
      ]
    ]
  }
}

1 feature, Polygon, 9 positions, 52,988 m²

All four return the same eight-cornered polygon covering 52,988 m². The two parcels are 28,136 m² and 28,135 m² on their own, so the union is 3,282 m² smaller than the two added together: the overlap is counted once.

Note the GDAL example needs two commands. ogr2ogr reads one dataset at a time, so the two parcels go into one layer before the union can see both. Use the GeoJSON union tool to merge any number of polygon features.

Intersection

Intersection returns the geometry common to both inputs, the area where they overlap. If the inputs do not overlap, the result is empty. Intersection quantifies how much of one area falls within another.

Real Scenario: Overlap Between Flood Zone and Residential Area

An emergency management agency needs to find the exact overlap between the FEMA 100-year flood zone and a residential zoning district. That polygon is the area both in the flood zone and zoned residential, the highest-priority area for flood insurance outreach.

Input 1: FEMA 100-year flood zone (a.geojson):

json{
  "type": "Feature",
  "properties": {
    "zone": "AE",
    "bfe_ft": 14
  },
  "geometry": {
    "type": "Polygon",
    "coordinates": [
      [
        [
          -122.482,
          37.818
        ],
        [
          -122.475,
          37.818
        ],
        [
          -122.475,
          37.822
        ],
        [
          -122.482,
          37.822
        ],
        [
          -122.482,
          37.818
        ]
      ]
    ]
  }
}

Input 2: Residential zoning district (b.geojson):

json{
  "type": "Feature",
  "properties": {
    "zone_code": "R-2",
    "units_per_acre": 8
  },
  "geometry": {
    "type": "Polygon",
    "coordinates": [
      [
        [
          -122.48,
          37.817
        ],
        [
          -122.473,
          37.817
        ],
        [
          -122.473,
          37.821
        ],
        [
          -122.48,
          37.821
        ],
        [
          -122.48,
          37.817
        ]
      ]
    ]
  }
}
python · Shapely 2.ximport json
from shapely.geometry import shape, mapping

flood = json.load(open("a.geojson"))
zoning = json.load(open("b.geojson"))

overlap = shape(flood["geometry"]).intersection(shape(zoning["geometry"]))

# An intersection can come back empty. Check before using it.
if overlap.is_empty:
    raise SystemExit("the two areas do not overlap")

print(json.dumps({"type": "Feature", "properties": {}, "geometry": mapping(overlap)}))

What it printed

{
  "type": "Feature",
  "properties": {},
  "geometry": {
    "type": "Polygon",
    "coordinates": [
      [
        [
          -122.475,
          37.818
        ],
        [
          -122.48,
          37.818
        ],
        [
          -122.48,
          37.821
        ],
        [
          -122.475,
          37.821
        ],
        [
          -122.475,
          37.818
        ]
      ]
    ]
  }
}

1 feature, Polygon, 5 positions, 146,590 m²

All four return the same rectangle, 146,590 m² or 14.7 ha. That number is the point of the operation: it is the area that is simultaneously in the flood zone and zoned residential, and it is what the report is written around.

An intersection can be empty, and every one of these examples has to handle that. Shapely returns an empty geometry, Turf returns null, PostGIS returns an empty GEOMETRYCOLLECTION, and ogr2ogr writes a file with no features. Two of the examples above check for it explicitly. Code that assumes an overlap exists fails on the first pair of polygons that do not touch.

Difference

Difference subtracts one polygon from another, returning the area of the first that the second does not cover. The operation is directional: A minus B differs from B minus A. If the second polygon lies entirely outside the first, the output is the first polygon unchanged.

Real Scenario: Subtract a Park From a City District

A zoning analyst needs the net developable area in a district after removing protected parkland. The inputs are the district boundary (A) and a park polygon (B) that sits wholly inside it. The output is the district with a hole in it.

Input A: City district (a.geojson):

json{
  "type": "Feature",
  "properties": {
    "district": "Midtown South"
  },
  "geometry": {
    "type": "Polygon",
    "coordinates": [
      [
        [
          -73.988,
          40.747
        ],
        [
          -73.982,
          40.747
        ],
        [
          -73.982,
          40.751
        ],
        [
          -73.988,
          40.751
        ],
        [
          -73.988,
          40.747
        ]
      ]
    ]
  }
}

Input B: Park to subtract (b.geojson):

json{
  "type": "Feature",
  "properties": {
    "name": "Bryant Park Annex",
    "protected": true
  },
  "geometry": {
    "type": "Polygon",
    "coordinates": [
      [
        [
          -73.9865,
          40.7482
        ],
        [
          -73.9845,
          40.7482
        ],
        [
          -73.9845,
          40.7498
        ],
        [
          -73.9865,
          40.7498
        ],
        [
          -73.9865,
          40.7482
        ]
      ]
    ]
  }
}
python · Shapely 2.ximport json
from shapely.geometry import shape, mapping

district = json.load(open("a.geojson"))
park = json.load(open("b.geojson"))

net = shape(district["geometry"]).difference(shape(park["geometry"]))

print(json.dumps({
    "type": "Feature",
    "properties": {"district": district["properties"]["district"]},
    "geometry": mapping(net),
}))

What it printed

{
  "type": "Feature",
  "properties": {
    "district": "Midtown South"
  },
  "geometry": {
    "type": "Polygon",
    "coordinates": [
      [
        [
          -73.988,
          40.751
        ],
        [
          -73.982,
          40.751
        ],
        [
          -73.982,
          40.747
        ],
        [
          -73.988,
          40.747
        ],
        [
          -73.988,
          40.751
        ]
      ],
      [
        [
          -73.9845,
          40.7482
        ],
        [
          -73.9845,
          40.7498
        ],
        [
          -73.9865,
          40.7498
        ],
        [
          -73.9865,
          40.7482
        ],
        [
          -73.9845,
          40.7482
        ]
      ]
    ]
  }
}

1 feature, Polygon, 10 positions, 195,073 m²

All four return a polygon with two rings. The first is the district outline, the second is the hole where the park was. The district covers 22.5 ha and the park 3.0 ha, leaving 19.5 ha developable.

Look at the ring direction in the output. RFC 7946 requires an exterior ring to run counterclockwise and a hole to run clockwise. The input district ring above is counterclockwise; the exterior ring that comes back is clockwise. Three of these four libraries are built on GEOS, which does not track the GeoJSON winding convention, so the result is valid geometry with the wrong winding for the specification, and almost no validator checks winding, so nothing downstream reports it.

Dissolve

Dissolve merges polygons that share a common attribute value, combining all features with the same property into one geometry per value. Union merges everything regardless of attributes; dissolve groups by a field first, producing one output feature per unique value.

Real Scenario: Combine States Into Regions

To map regional data you dissolve state polygons grouped by a region property. The input below has four squares in two regions. Two of the three West squares share an edge and merge into one shape; the third sits apart. That is the case worth testing, because the result has to hold a region that is in two separate pieces.

Input: Four features in two regions (a.geojson):

json{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "properties": {
        "name": "West A",
        "region": "West"
      },
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [
              -120,
              40
            ],
            [
              -119,
              40
            ],
            [
              -119,
              41
            ],
            [
              -120,
              41
            ],
            [
              -120,
              40
            ]
          ]
        ]
      }
    },
    {
      "type": "Feature",
      "properties": {
        "name": "West B",
        "region": "West"
      },
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [
              -119,
              40
            ],
            [
              -118,
              40
            ],
            [
              -118,
              41
            ],
            [
              -119,
              41
            ],
            [
              -119,
              40
            ]
          ]
        ]
      }
    },
    {
      "type": "Feature",
      "properties": {
        "name": "West C",
        "region": "West"
      },
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [
              -116,
              40
            ],
            [
              -115,
              40
            ],
            [
              -115,
              41
            ],
            [
              -116,
              41
            ],
            [
              -116,
              40
            ]
          ]
        ]
      }
    },
    {
      "type": "Feature",
      "properties": {
        "name": "East A",
        "region": "East"
      },
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [
              -75,
              40
            ],
            [
              -74,
              40
            ],
            [
              -74,
              41
            ],
            [
              -75,
              41
            ],
            [
              -75,
              40
            ]
          ]
        ]
      }
    }
  ]
}
python · Shapely 2.ximport json
import geopandas as gpd

states = gpd.read_file("a.geojson")

# dissolve() groups by the column and unions each group. aggfunc controls
# what happens to the other columns.
regions = states.dissolve(by="region", aggfunc={"name": "count"})
regions = regions.rename(columns={"name": "feature_count"}).reset_index()

print(regions.to_json())

What it printed

{
  "type": "FeatureCollection",
  "features": [
    {
      "id": "0",
      "type": "Feature",
      "properties": {
        "region": "East",
        "feature_count": 1
      },
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [
              -75,
              40
            ],
            [
              -74,
              40
            ],
            [
              -74,
              41
            ],
            [
              -75,
              41
            ],
            [
              -75,
              40
            ]
          ]
        ]
      }
    },
    {
      "id": "1",
      "type": "Feature",
      "properties": {
        "region": "West",
        "feature_count": 3
      },
      "geometry": {
        "type": "MultiPolygon",
        "coordinates": [
          [
            [
              [
                -120,
                40
              ],
              [
                -120,
                41
              ],
              [
                -119,
                41
              ],
              [
                -118,
                41
              ],
              [
                -118,
                40
              ],
              [
                -119,
                40
              ],
              [
                -120,
                40
              ]
            ]
          ],
          [
            [
              [
                -115,
                41
              ],
              [
                -115,
                40
              ],
              [
                -116,
                40
              ],
              [
                -116,
                41
              ],
              [
                -115,
                41
              ]
            ]
          ]
        ]
      }
    }
  ]
}

2 features, MultiPolygon and Polygon, 17 positions, 37651.155 km²

Three of the four return two features and Turf returns three. Shapely, PostGIS and GDAL each give one feature per region, and the West region comes back as a MultiPolygon with two parts: the merged pair, and the square that touches nothing. Turf keeps every merged piece as its own feature, so the West region arrives as two separate Polygons that both carry region: West.

The ground covered is the same either way. The structure is not, and code that expects one feature per group gets three where it expected two. Group Turf's output yourself if you need one feature per value.

Combining Operations: A Real Workflow

Individual operations become useful when chained. A common workflow in urban planning, environmental analysis, and logistics combines three or more to derive a result no single operation could produce.

Workflow: Residential Land Within Walking Distance of a Transit Stop

Goal: how much residential-zoned land in a district falls within 800m, about a 10-minute walk, of a transit station entrance.

  1. Buffer the station point by 800m, using one of the four approaches above. Remember the projection step, and remember the four will not agree to the metre.
  2. Intersect that buffer with the district boundary, which trims the walking-distance circle to the part inside the district.
  3. Clip the residential parcels to the result of step 2, keeping only the parts inside both the walk zone and the district.
  4. Sum the area of the clipped parcels. On the WGS 84 ellipsoid, not in square degrees: ST_Area(geog) in PostGIS, turf.area() in JavaScript, or pyproj.Geod.geometry_area_perimeter() in Python.

The answer is one number, and the buffer spread above is the reason to quote it with a margin. A 0.12% difference at step 1 carries through every later step, so two teams running this workflow in different tools will not get the same hectare count. The pattern recurs across transit planning, retail site selection, emergency response and environmental impact work.

Versions used: shapely 2.1.2, geopandas 1.1.4, gdal 3.8.4, postgis 3.4.2, turf 7.3.4. The results are published under CC BY 4.0 on the datasets page.