GeoPandas Basics: Working with Spatial Data in Python
On this page
Problem statement
A first session with GeoPandas is the same five moves whatever the data: read a file, look at what is in it, check its coordinate system, keep the rows you need, and write the result. Each move is one line. What makes the session go right or wrong is what the file turns out to hold — and that only shows if you look at the output of each line rather than past it.
This page makes those five moves on one small real file, OpenStreetMap's land use for Rutland, England's smallest historic county, and stops at each to say what it found. Every step links to the guide that covers it in depth.
Quick answer
Read, inspect, check the CRS, filter, save — five lines on a GeoDataFrame, a pandas table with a geometry column (what a GeoDataFrame is):
import geopandas as gpd
landuse = gpd.read_file("data/rutland/gis_osm_landuse_a_free_1.shp")
print(len(landuse), "polygons in", landuse.crs)
residential = landuse[landuse["fclass"] == "residential"]
residential.to_file("output/residential.geojson", driver="GeoJSON")
print(len(residential), "residential polygons written")2405 polygons in EPSG:4326
159 residential polygons written
The file is OpenStreetMap land use for Rutland, England's smallest historic county, from the Geofabrik shapefile extract; the section The data used on this page downloads it. Every output on this page was printed by running the code.
The vector workflow pipeline
Read, inspect, check the CRS, filter, plot, save
Read a shapefile or GeoJSON into a GeoDataFrame
gpd.read_file() reads both formats used here (installing it: setting up a Python GIS environment).
Read a shapefile
A shapefile is several files side by side (.shp, .shx, .dbf, .prj, often .cpg); you point read_file() at the .shp and it finds the rest.
import geopandas as gpd
landuse = gpd.read_file("data/rutland/gis_osm_landuse_a_free_1.shp")
print(landuse.head()) osm_id code fclass name geometry
0 3554658 7210 nature_reserve Lagoon 3 POLYGON ((-0.69395 52.6626, -0.69361...
1 3554659 7210 nature_reserve Lagoon 2 POLYGON ((-0.70015 52.65614, -0.6991...
2 3700148 7202 park NaN POLYGON ((-0.74126 52.66536, -0.7410...
3 3700814 7202 park Royce Recreation Ground POLYGON ((-0.74212 52.67041, -0.7405...
4 3700816 7202 park NaN POLYGON ((-0.74291 52.66793, -0.7424...
Read a GeoJSON file
The second layer is Rutland's 15 electoral wards from the ONS boundary service, saved as GeoJSON.
wards = gpd.read_file("data/rutland_wards.geojson")
print(wards.head()) WD24CD WD24NM geometry
0 E05012632 Barleythorpe POLYGON ((-0.72525 52.68542, -0.7243...
1 E05012633 Braunston & Martinsthorpe POLYGON ((-0.75976 52.65781, -0.7594...
2 E05012634 Cottesmore POLYGON ((-0.67132 52.75641, -0.6707...
3 E05012635 Exton POLYGON ((-0.70772 52.72347, -0.709 ...
4 E05012636 Greetham POLYGON ((-0.60859 52.75971, -0.6054...
Each call returns a GeoDataFrame; the details of each format are in reading a shapefile and reading GeoJSON.
Inspect the GeoDataFrame structure
Look before using anything:
print(landuse.columns.tolist())
print(landuse.dtypes)
print("geometry column:", landuse.geometry.name)
print(len(landuse), "rows")['osm_id', 'code', 'fclass', 'name', 'geometry']
osm_id str
code int32
fclass str
name str
geometry geometry
dtype: object
geometry column: geometry
2405 rows
Five columns and 2,405 rows. Geofabrik's code is a numeric version of the fclass text, and one column is mostly empty:
print(landuse["name"].isna().sum(), "of", len(landuse), "polygons have no name")2165 of 2405 polygons have no name
A filter on name would silently skip 91% of the county's land use. Exploring an unfamiliar dataset turns this look into a routine.
Check the coordinate reference system
Both layers report their CRS:
print(landuse.crs)
print(wards.crs)EPSG:4326
EPSG:4326
Both layers are in longitude and latitude (WGS 84): Geofabrik ships every extract that way, and the ward download asked the ONS service for it. That is fine for reading and plotting, but the units are degrees, so an area computed directly is a number of square degrees. Reproject to a metric CRS — British National Grid here — before measuring:
print(f"{landuse.to_crs(27700).area.sum() / 1e6:,.1f} km² of mapped land use")268.9 km² of mapped land use
Why degrees cannot be measured is projected vs geographic CRS.
Filter spatial records
A common task is filtering features based on an attribute field. Look at what values the field holds first:
print(landuse["fclass"].value_counts().head(8))fclass
forest 1247
grass 323
farmland 241
residential 159
scrub 93
farmyard 90
park 40
meadow 37
Name: count, dtype: int64
Example: keep only residential land
residential = landuse[landuse["fclass"] == "residential"]
print(len(residential), "residential polygons")
print(residential["name"].dropna().head().tolist())159 residential polygons
['Wing', 'Cottesmore Hall Redevelopment', "Permanent Traveller's Camp", 'Manor Farm', 'Lonsdale Park']
You can also remove rows with missing or empty geometry:
with_geometry = landuse[landuse.geometry.notna() & ~landuse.geometry.is_empty]
print(len(landuse) - len(with_geometry), "rows dropped")0 rows dropped
None here; what null, empty and invalid geometry are, and why they differ, is in null, empty and invalid geometry.
Check geometry types
print(landuse.geom_type.value_counts())
print(wards.geom_type.value_counts())Polygon 2401
MultiPolygon 4
Name: count, dtype: int64
Polygon 15
Name: count, dtype: int64
Four land-use areas are drawn in more than one piece and come back as MultiPolygon; a shapefile does not distinguish the two, so a polygon layer is often a mix. That is normal. Points or lines in a polygon layer are not.
Plot the data for a quick visual check
A plot is the fastest check that the layer is where it should be:
residential.plot()For a slightly clearer result, draw the wards underneath so you can see whether the filtered features land where you expect:
ax = wards.plot(figsize=(8, 8), color="white", edgecolor="grey")
residential.plot(ax=ax, color="tab:red")
ax.set_title("Residential land use in Rutland")A layer in the wrong CRS shows up at once: it lands far away from the other one, or the axes read in hundreds of thousands instead of degrees. Making the map presentable is plotting maps with GeoPandas and Matplotlib.
Save the result to a new file
Save to GeoJSON
residential.to_file("output/residential.geojson", driver="GeoJSON")Save to a shapefile or GeoPackage
residential.to_file("output/residential.shp")
residential.to_file("output/rutland.gpkg", layer="residential")
import os
shp = sorted(f for f in os.listdir("output") if f.startswith("residential.") and not f.endswith((".geojson", ".png")))
print(shp, f"{sum(os.path.getsize('output/' + f) for f in shp) / 1e3:.1f} kB")
print(f"residential.geojson {os.path.getsize('output/residential.geojson') / 1e3:.1f} kB")
print(f"rutland.gpkg {os.path.getsize('output/rutland.gpkg') / 1e3:.1f} kB")['residential.cpg', 'residential.dbf', 'residential.prj', 'residential.shp', 'residential.shx'] 138.8 kB
residential.geojson 181.6 kB
rutland.gpkg 225.3 kB
One shapefile write is five files on disk. Sizes differ too: the .dbf stores every text field at a fixed width — 80 characters for each of the three here — so it is 40 kB for 159 rows, 250 bytes a row however short the values, while GeoJSON writes each coordinate out as text and a GeoPackage is an SQLite database whose tables and spatial index cost about 98 kB before any data — a one-feature GeoPackage of this layer was 98.3 kB. None of this matters at this size, but it is worth knowing which way each format grows; GIS file formats compared goes further.
Rutland land use summarised, filtered and exported
The data used on this page
Two real layers, both for Rutland:
- Land use —
gis_osm_landuse_a_free_1.shpfrom Geofabrik's daily OpenStreetMap extract (© OpenStreetMap contributors, ODbL); counts drift a little as the map is edited. - Wards — Rutland's 15 wards from the ONS December 2024 boundaries (OGL), as GeoJSON.
import io
import os
import zipfile
import geopandas as gpd
import pandas as pd
import requests
pd.set_option("display.width", 120) # print whole rows instead of eliding columns
pd.set_option("display.max_columns", None)
pd.set_option("display.max_colwidth", 40)
os.makedirs("data/rutland", exist_ok=True)
os.makedirs("output", exist_ok=True)
if not os.path.exists("data/rutland/gis_osm_landuse_a_free_1.shp"):
extract = requests.get("https://download.geofabrik.de/europe/united-kingdom/england/"
"rutland-260101-free.shp.zip", timeout=300)
extract.raise_for_status()
zipfile.ZipFile(io.BytesIO(extract.content)).extractall("data/rutland")
if not os.path.exists("data/rutland_wards.geojson"):
wards = requests.get("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
"Wards_December_2024_Boundaries_UK_BGC/FeatureServer/0/query",
params={"where": "LAD24NM='Rutland'", "outFields": "WD24CD,WD24NM",
"outSR": 4326, "f": "geojson"}, timeout=120)
wards.raise_for_status()
open("data/rutland_wards.geojson", "wb").write(wards.content)
print(len([f for f in os.listdir("data/rutland") if f.endswith(".shp")]), "shapefiles in the extract")18 shapefiles in the extract
The extract holds a shapefile per theme — roads, buildings, water, places and so on. This page reads only the land-use one.
Example 1: Load a shapefile and summarise it
import geopandas as gpd
landuse = gpd.read_file("data/rutland/gis_osm_landuse_a_free_1.shp")
print(f"{len(landuse):,} features, CRS {landuse.crs}, columns {landuse.columns.tolist()}")
print("bounds:", landuse.total_bounds.round(3))2,405 features, CRS EPSG:4326, columns ['osm_id', 'code', 'fclass', 'name', 'geometry']
bounds: [-0.855 52.53 -0.395 52.763]
total_bounds is [minx, miny, maxx, maxy] — here, longitudes and latitudes spanning the county, about 0.46° by 0.23°. It is the quickest check that a layer is where you think it is.
Example 2: Load a GeoJSON file and inspect columns
import geopandas as gpd
wards = gpd.read_file("data/rutland_wards.geojson")
print(wards.columns.tolist())
print(wards.geometry.name)
print(wards[["WD24CD", "WD24NM"]].sort_values("WD24NM").head())['WD24CD', 'WD24NM', 'geometry']
geometry
WD24CD WD24NM
0 E05012632 Barleythorpe
1 E05012633 Braunston & Martinsthorpe
2 E05012634 Cottesmore
3 E05012635 Exton
4 E05012636 Greetham
The geometry column is called geometry whatever the file itself called it — a GeoPackage column named geom also reads back as geometry — and the ward codes are the ONS E05… codes, which are what you would join population or election tables on.
Example 3: Check CRS and geometry types
import geopandas as gpd
wards = gpd.read_file("data/rutland_wards.geojson")
print("CRS:", wards.crs, "| geographic:", wards.crs.is_geographic)
print(wards.geom_type.value_counts())CRS: EPSG:4326 | geographic: True
Polygon 15
Name: count, dtype: int64
crs.is_geographic is True for a longitude/latitude CRS — the signal to reproject before you measure lengths or areas.
Example 4: Filter features by attribute
import geopandas as gpd
landuse = gpd.read_file("data/rutland/gis_osm_landuse_a_free_1.shp")
commercial = landuse[landuse["fclass"].isin(["commercial", "retail", "industrial"])]
print(commercial["fclass"].value_counts())
print(commercial.loc[commercial["name"].notna(), ["fclass", "name"]].head())fclass
industrial 30
commercial 8
retail 6
Name: count, dtype: int64
fclass name
19 commercial Oakham Enterprise Park
205 industrial Wing Water Treatment Works
461 industrial Ketton Cement Works
749 industrial Tixover
833 industrial Glebe Works
isin() keeps any of several values in one step; chaining == comparisons with | does the same thing but is easier to get wrong.
Example 5: Plot and export filtered data
import geopandas as gpd
landuse = gpd.read_file("data/rutland/gis_osm_landuse_a_free_1.shp")
wards = gpd.read_file("data/rutland_wards.geojson")
residential = landuse[(landuse["fclass"] == "residential") & landuse.geometry.notna()]
ax = wards.plot(figsize=(8, 8), color="white", edgecolor="grey")
residential.plot(ax=ax, color="tab:red")
ax.set_title("Residential land use in Rutland")
ax.figure.savefig("output/residential.png", dpi=100)
residential.to_file("output/residential.geojson", driver="GeoJSON")
print(len(gpd.read_file("output/residential.geojson")), "features read back from the GeoJSON")159 features read back from the GeoJSON
Reading the output back is the cheapest test that the write did what you meant.
What one small file had to say
Every step on this page is one line, and every one of them succeeded. What made them worth doing was reading their output. The column listing showed that 2,165 of 2,405 polygons have no name — invisible to anyone who went straight to filtering. The CRS line said degrees, so the area had to wait for a reprojection before it meant anything (268.9 km² of mapped land use). The geometry types showed four multipolygons in a "polygon" layer. The write showed that one shapefile is five files and that its text fields take the same space however short their values.
None of those raised an error, and none needed more than one line to see.
Invalid geometry, missing labels and big files
Missing or invalid geometry
Some rows may have null geometry values or invalid shapes. Null geometry is dropped with a filter; validity is a separate check:
print(landuse.geometry.isna().sum(), "null geometries")
print(landuse.geometry.is_valid.value_counts())0 null geometries
True 2405
Name: count, dtype: int64
This extract is clean. Data digitised by hand or exported from a CAD system often is not, and invalid geometry makes overlays and dissolves fail or return wrong areas. How to fix invalid geometries covers the repair.
CRS missing, field names cut short
Copy a shapefile without its .prj and it reads with no CRS; label it only from what the source documents (cannot transform naive geometries). Write a column called area_hectares to a shapefile and it comes back as area_hecta, because the format allows ten characters (shapefile column names truncated).
Large files can be slow
Large parcel or boundary datasets may load slowly and be expensive to plot. You do not have to read a whole file to see what is in it — pyogrio.read_info() reports the count, CRS and fields without loading the features, and rows= reads just the first few:
import pyogrio
info = pyogrio.read_info("data/rutland/gis_osm_buildings_a_free_1.shp")
print(info["features"], "features;", info["crs"], ";", list(info["fields"]))
print(gpd.read_file("data/rutland/gis_osm_buildings_a_free_1.shp", rows=3)[["osm_id", "fclass"]])5513 features; EPSG:4326 ; ['osm_id', 'code', 'fclass', 'name', 'type']
osm_id fclass
0 3700606 building
1 3700607 building
2 59833301 building
read_file() also takes columns= to skip fields and bbox= or mask= to read only the features in an area.
Internal links
- Python for GIS: What It Is and When to Use It
- Coordinate Reference Systems (CRS) Explained for Python GIS
- Overlay Operations in GeoPandas: Union, Intersection, Difference Explained
- Geoprivacy explained: why coordinates are personal data — when a geometry column is personal data
FAQ
What is GeoPandas used for?
GeoPandas is used for working with vector spatial data in Python. Common tasks include reading shapefiles and GeoJSON, filtering features, checking CRS, plotting data, and exporting results.
What is the difference between pandas and GeoPandas?
pandas handles standard tabular data. GeoPandas adds a geometry column so the table can store and work with spatial features such as points, lines, and polygons.
How do I select specific columns from a GeoDataFrame?
Index it like a pandas DataFrame, but keep the geometry column if you want to stay spatial: gdf[["osm_id", "fclass", "geometry"]]. Dropping the geometry column returns a plain DataFrame rather than a GeoDataFrame.
How do I count features by attribute value in GeoPandas?
Use value_counts() on the column, for example gdf["fclass"].value_counts(). This works exactly as in pandas and is a quick way to summarize categories before filtering or plotting.
Can I create a GeoDataFrame from a CSV with latitude and longitude columns?
Yes. Build point geometries from the coordinate columns and set the CRS explicitly:
gdf = gpd.GeoDataFrame(
df, geometry=gpd.points_from_xy(df["lon"], df["lat"]), crs="EPSG:4326"
)