Shapely Basics: Working with Geometry Objects in Python

On this page

Problem statement

Shapely is the geometry underneath GeoPandas: points from GPS coordinates, lines for routes, polygons for parcels, and tests of whether one is inside another. It gives you objects you can create, inspect, compare, and transform in memory. The examples below build real ones: the Black Bull Inn at Market Overton in Rutland as a point, Teigh Road outside it as a line, and the pub's Land Registry title parcel as a polygon — all in British National Grid metres, so lengths and areas mean something.

This page covers:

creating common Shapely geometry objects
reading basic geometry properties
checking spatial relationships
running simple geometry operations
converting geometry to readable text

This page does not cover:

reading shapefiles or GeoJSON with GeoPandas
CRS management in detail
advanced topology rules or performance tuning

Quick answer

The three core Shapely geometry types: Point, LineString, and Polygon.

Shapely lets you work with geometry objects such as Point, LineString, and Polygon in pure Python.

Basic workflow:

import Shapely geometry classes
create objects from coordinates
inspect properties like area, length, bounds, and centroid
test relationships like intersects() and contains()
apply simple operations like buffer()

from shapely.geometry import LineString, Point, Polygon

pub = Point(488707.4, 316459.6)                 # the Black Bull Inn, EPSG:27700
road = LineString(TEIGH_ROAD)                   # Teigh Road's vertices from OpenStreetMap
parcel = Polygon(PARCEL_RING)                   # Land Registry title parcel 31894662

print(pub.geom_type, "|", f"{road.length:.1f} m of road |", f"{parcel.area:.1f} m² of parcel")
print("parcel contains pub:", parcel.contains(pub))

Point | 739.1 m of road | 292.0 m² of parcel
parcel contains pub: True

Points, lines and polygons, their properties and relations

Each type is built from the one above it, which is why the constructors nest.

Import Shapely geometry classes

For basic geometry work, import the main geometry types from shapely.geometry:

from shapely.geometry import LineString, MultiLineString, MultiPoint, MultiPolygon, Point, Polygon

The most common classes are:

Point for a single location
LineString for a connected path
Polygon for an area
multipart types such as MultiPoint and MultiPolygon when one feature has multiple parts

Create a Point geometry

A Point is built from x and y coordinates:

site = Point(488707.4, 316459.6)
print(site, "|", site.geom_type, "|", site.x, site.y)

POINT (488707.4 316459.6) | Point | 488707.4 316459.6

Typical GIS uses are site locations, GPS observations, address coordinates and sampling points. This is the standard pattern when you need to create geometry from raw coordinates.

Create a LineString geometry

A LineString is an ordered sequence of coordinate pairs:

print(len(TEIGH_ROAD), "vertices, the first three:", TEIGH_ROAD[:3])
road_segment = LineString(TEIGH_ROAD)
print(road_segment.geom_type, f"{road_segment.length:.1f} m")
print("reversed is equal:", road_segment.equals(LineString(TEIGH_ROAD[::-1])),
      "| same coordinates in the same order:", road_segment.equals_exact(LineString(TEIGH_ROAD[::-1]), 0))

14 vertices, the first three: [(488716.74, 316474.89), (488679.48, 316456.53), (488634.67, 316423.35)]
LineString 739.1 m
reversed is equal: True | same coordinates in the same order: False

Order matters because the line follows the coordinates in the sequence you provide — equals treats a reversed line as the same set of points, equals_exact does not. Typical uses are roads, rivers, walking paths and route segments.

Create a Polygon geometry

A Polygon represents an area, created from a ring of coordinates:

print(len(PARCEL_RING), "coordinates; first == last:", PARCEL_RING[0] == PARCEL_RING[-1])
study_area = Polygon(PARCEL_RING)
unclosed = Polygon(PARCEL_RING[:-1])
print(study_area.geom_type, "| closed for you:", unclosed.equals(study_area))

17 coordinates; first == last: True
Polygon | closed for you: True

Shapely closes the ring for you if the last coordinate is missing, but repeating the first coordinate at the end makes the definition clearer. Typical uses are land parcels, administrative boundaries, service zones and site extents.

Inspect geometry properties

Once you create a geometry, you can inspect its key properties:

print(study_area.geom_type)                                   # Polygon
print(tuple(round(v, 1) for v in study_area.bounds))          # (minx, miny, maxx, maxy)
print(study_area.centroid)                                    # centre point
print(round(study_area.area, 1), "m²")                        # area for polygons
print(round(study_area.length, 1), "m perimeter")             # length of a polygon is its perimeter
print(round(road_segment.length, 1), "m of road")             # line length

Polygon
(488690.6, 316448.0, 488717.9, 316469.7)
POINT (488704.29475896887 316458.6838492458)
292.0 m²
72.6 m perimeter
739.1 m of road

Useful properties in GIS work:

geom_type: confirms the geometry type
bounds: helpful for extent checks and filtering
centroid: useful for labels or quick reference points
area: used for polygon measurement
length: line length, or a polygon's perimeter
point x and y: direct coordinate access

The properties above, from real coordinates.

Check spatial relationships between geometries

Shapely supports the common spatial relationship checks:

print("parcel contains the pub:", study_area.contains(site))
print("pub within the parcel:  ", site.within(study_area))
print("road intersects parcel: ", road_segment.intersects(study_area))
print("road touches parcel:    ", road_segment.touches(study_area))
print(f"road to parcel: {road_segment.distance(study_area):.2f} m")

parcel contains the pub: True
pub within the parcel:   True
road intersects parcel:  False
road touches parcel:     False
road to parcel: 2.52 m

These methods are useful for point-in-polygon tests, overlap checks, boundary contact checks and filtering features by location. The road runs past the parcel without touching it: Land Registry parcels usually stop at the property boundary and OpenStreetMap draws the road along its centre line, so a gap of a few metres is normal, not an error.

Run simple geometry operations

Shapely geometry operations return new geometry objects:

buffered_site = site.buffer(100)
print(buffered_site.geom_type, f"{buffered_site.area:,.0f} m² (a true circle would be {3.14159265 * 100**2:,.0f})")

intersection_geom = road_segment.intersection(buffered_site)
print(intersection_geom.geom_type, f"{intersection_geom.length:.1f} m of road within 100 m")

print(f"distance from the pub to the parcel boundary: {site.distance(study_area.boundary):.2f} m")

combined = study_area.union(buffered_site)
print(combined.geom_type, f"{combined.area:,.0f} m²")

Polygon 31,365 m² (a true circle would be 31,416)
LineString 116.2 m of road within 100 m
distance from the pub to the parcel boundary: 5.47 m
Polygon 31,365 m²

buffer() approximates a circle with straight segments — 16 per quarter circle by default — so its area is slightly under πr². Common operations:

buffer(distance) creates an area around a geometry
intersection(other) returns the shared part
union(other) merges geometries
distance(other) calculates distance between geometries

Convert geometry to readable formats

Readable output helps with debugging and exporting:

print(site.wkt)
print(road_segment.wkt[:70], "…")
print(site.__geo_interface__)

POINT (488707.4 316459.6)
LINESTRING (488716.74 316474.89, 488679.48 316456.53, 488634.67 316423 …
{'type': 'Point', 'coordinates': (488707.4, 316459.6)}

WKT (POINT (488707.4 316459.6)) is useful for logs, debugging, database inserts and quick inspection; __geo_interface__ is a GeoJSON-like mapping.

Multipart parts, merges, WKT and a malformed polygon

The data used on this page

The Black Bull Inn and Teigh Road from Geofabrik's OpenStreetMap extract of 1 January 2026 (© OpenStreetMap contributors, ODbL), and title parcel 31894662 from HM Land Registry's INSPIRE parcels for Rutland (© Crown copyright and database rights; OGL). The setup extracts their coordinates as plain Python lists, which is the form coordinates arrive in from most sources.

from pathlib import Path

import geopandas as gpd
import requests
import shapely

Path("data").mkdir(exist_ok=True)
for target, url in (("data/rutland.shp.zip", "https://download.geofabrik.de/europe/united-kingdom/england/rutland-260101-free.shp.zip"),
                    ("data/Rutland.zip", "https://use-land-property-data.service.gov.uk/datasets/inspire/download/Rutland.zip")):
    if not Path(target).exists():
        Path(target).write_bytes(requests.get(url, timeout=300).content)

parcel_row = gpd.read_file("/vsizip/data/Rutland.zip/Land_Registry_Cadastral_Parcels.gml",
                           where="INSPIREID = 31894662", DOWNLOAD_SCHEMA="NO")
roads = gpd.read_file("/vsizip/data/rutland.shp.zip/gis_osm_roads_free_1.shp", where="name = 'Teigh Road'").to_crs(27700)
nearest = roads.distance(parcel_row.geometry.iloc[0]).idxmin()

PARCEL_RING = [(round(x, 2), round(y, 2)) for x, y in parcel_row.geometry.iloc[0].exterior.coords]
TEIGH_ROAD = [(round(x, 2), round(y, 2)) for x, y in roads.geometry.loc[nearest].coords]
print("Shapely", shapely.__version__, "| parcel ring:", len(PARCEL_RING), "coordinates | road:", len(TEIGH_ROAD), "vertices")

Shapely 2.1.2 | parcel ring: 17 coordinates | road: 14 vertices

Example 1: Multipart geometries and their parts

pubs = gpd.read_file("/vsizip/data/rutland.shp.zip/gis_osm_pois_free_1.shp", where="fclass = 'pub'").to_crs(27700)
all_pubs = MultiPoint(list(pubs.geometry))
print(all_pubs.geom_type, len(all_pubs.geoms), "parts")
print("first part:", all_pubs.geoms[0])
try:
    all_pubs[0]
except TypeError as exc:
    print("indexing the multipart directly:", exc)

MultiPoint 33 parts
first part: POINT (490479.934860832 298301.25909281796)
indexing the multipart directly: 'MultiPoint' object is not subscriptable

Iterate or index through .geoms; Shapely 2 removed direct indexing of multipart geometries.

Example 2: Merge many geometries at once

from shapely import union_all

buffers = [p.buffer(500) for p in pubs.geometry]
merged = union_all(buffers)
print(len(buffers), "500 m buffers ->", merged.geom_type, "with", len(merged.geoms), "parts,",
      f"{merged.area / 1e6:.1f} km² covered")

33 500 m buffers -> MultiPolygon with 23 parts, 20.9 km² covered

union_all dissolves a list of geometries in one call; overlapping buffers in the same village merge into one part.

Example 3: Round-trip through WKT

from shapely import from_wkt

text = study_area.wkt
back = from_wkt(text)
print(len(text), "characters of WKT | equal after round trip:", back.equals(study_area))

364 characters of WKT | equal after round trip: True

Example 4: A malformed Polygon call

try:
    Polygon((0, 0), (10, 0), (10, 10))
except Exception as exc:
    print(type(exc).__name__, "—", str(exc)[:90])
print(Polygon([(0, 0), (10, 0), (10, 10)]).area)

TypeError — Polygon.__new__() takes from 1 to 3 positional arguments but 4 were given
50.0

A Polygon's first argument is the whole exterior ring; the only other argument is a list of holes. Three points passed as separate arguments are one argument too many, so the call fails before Shapely looks at a coordinate. Pass one list — Shapely closes the ring and the 10 × 10 right triangle has an area of 50.

Shapely is geometry, without files or attributes

Shapely works with geometry objects in memory. It does not read shapefiles by itself and it does not manage attribute tables like GeoPandas does. Its job is geometry creation and geometry logic.

Geometries are immutable, so the result of an operation must be assigned.

It helps to separate three things:

coordinate lists: raw tuples such as TEIGH_ROAD
geometry objects: Shapely objects such as LineString(TEIGH_ROAD)
geometry methods and properties: tools like .area, .buffer(), and .intersects()

In practice, many Python GIS workflows start here. You create or receive geometry objects, test them, measure them, or transform them, then later place them into a GeoDataFrame for file export, filtering, clipping, joins, or larger analysis.

No CRS, valid rings and coordinate units

Geometry objects do not store CRS by themselves

A Shapely geometry is just coordinates. It does not know whether those values are in EPSG:4326, UTM, or a local projected CRS. If you need CRS-aware workflows, use GeoPandas or manage CRS separately.

Polygon rings should be valid

Self-intersecting polygons or badly constructed rings can produce invalid geometry. That can break overlay or area calculations. You can do a basic check with geometry.is_valid.

print(study_area.is_valid)

True

If you run into problems, see How to Fix Invalid Geometry Errors in Shapely.

Area and length depend on coordinate units

If your coordinates are in longitude and latitude, area and length are not in meters. For meaningful measurements, use an appropriate projected CRS before calculating. For background, see Coordinate Reference Systems (CRS) Explained for Python GIS.

Operations return new geometries

Methods such as buffer() and intersection() do not change the original object in place.

new_geom = site.buffer(100)

Some results may be multipart geometries

An intersection or union can return MultiPolygon, MultiLineString, or an empty geometry depending on the inputs.

Internal links

For broader context, read Python for GIS: What It Is and When to Use It.

Related background pages:

Troubleshooting:

How to Fix Invalid Geometry Errors in Shapely

FAQ

What is Shapely used for in Python GIS?

Shapely is used to create and manipulate geometry objects in Python. Common tasks include building points, lines, and polygons, checking intersections, measuring area or distance, and running operations such as buffer or union.

What geometry types can I create with Shapely?

You can create Point, LineString, Polygon, MultiPoint, MultiLineString, and MultiPolygon, along with other geometry types such as geometry collections.

Does Shapely handle coordinate reference systems?

No. Shapely geometry objects do not store CRS information by themselves. If you need CRS management, use GeoPandas or another GIS structure that tracks CRS separately.

What is the difference between Shapely and GeoPandas?

Shapely handles geometry objects and geometry operations. GeoPandas builds on Shapely and adds tabular data handling, file I/O, CRS support, and dataframe-style spatial workflows.

Why is my polygon area incorrect in Shapely?

The most common reason is using geographic coordinates instead of a projected CRS. If your polygon coordinates are in degrees, the area result is not a real square-meter value. Invalid polygon geometry can also cause incorrect results.

How do I build many Shapely points from columns of x and y values?

For a single geometry use Point(x, y), but to vectorize a whole table use GeoPandas: gpd.points_from_xy(df["lon"], df["lat"]) returns a geometry array you can assign to a GeoDataFrame. This is much faster than looping and creating one Point at a time.

How do I merge a list of geometries into one in modern Shapely?

Use shapely.union_all(geometries), the Shapely 2 function for what shapely.ops.unary_union does; the older name still works, without a warning, in Shapely 2.1. It dissolves a list or array of geometries into a single combined geometry in one call.

How do I access the individual parts of a MultiPolygon or GeometryCollection?

Iterate over the .geoms attribute, for example for part in multi.geoms: print(part.area). Indexing a multipart geometry directly is not supported; you must go through .geoms.

How do I create a Shapely geometry from a WKT string?

Use shapely.from_wkt, for example from_wkt("POINT (12.5 41.9)"), which parses the text back into a geometry object. The reverse direction is the .wkt property on any geometry.

Why does `Polygon` give a coordinate or "not enough values" error when I pass it points?

A Polygon's first argument must be a sequence of coordinate tuples (the exterior ring), not separate arguments, so Polygon((0,0),(1,0),(1,1)) fails. Pass a single list instead: Polygon([(0,0),(1,0),(1,1),(0,0)]).