--- title: "Getting started with GeoPandas: projections, spatial joins and choropleth maps" description: "Learn GeoPandas with Natural Earth data: reproject to an equal-area CRS, map population density, join cities to countries and dissolve to continents." url: https://www.technicalgis.com/tutorials/geopandas-getting-started/ author: TechnicalGIS Team publisher: TechnicalGIS (https://www.technicalgis.com/about/) published: 2026-10-14 updated: 2026-10-14 language: en-GB --- # Getting started with GeoPandas: projections, spatial joins and choropleth maps > Learn GeoPandas with Natural Earth data: reproject to an equal-area CRS, map population density, join cities to countries and dissolve to continents. - Level: Beginner - Tested with: Python 3.8 - Source code: [Download scripts](https://www.technicalgis.com/code/geopandas-getting-started/) ## Key takeaways 1. A GeoDataFrame is a pandas DataFrame with a geometry column and a CRS, so vector GIS operations sit alongside ordinary pandas code. 2. Reproject before you measure: areas in EPSG:4326 are in square degrees, so use an equal-area projection such as World Mollweide for global areas. 3. sjoin() attaches attributes by location. Count the unmatched rows afterwards: 7% of Natural Earth cities, mostly coastal, missed a country until sjoin\_nearest() was used. 4. dissolve() is a spatial groupby, and plot(column=..., scheme=...) draws a classified choropleth in one call. ## Before you start - Python 3.8 or newer with pip and a virtual environment - Basic pandas: selecting columns, groupby - Around 10 MB of disk space for the Natural Earth downloads If you have used QGIS or ArcGIS Pro, you already know the core operations of vector GIS: reprojecting, calculating areas, joining layers by location and dissolving boundaries. [GeoPandas](https://geopandas.org/) gives you the same operations in Python. Because a `GeoDataFrame` is a pandas `DataFrame` with a geometry column, your analysis becomes a script you can re-run, review and schedule. In this tutorial we will use free [Natural Earth](https://www.naturalearthdata.com/) data to: 1. read country polygons and city points straight from a zip file; 2. reproject them to an **equal-area** coordinate reference system (CRS) so areas are correct; 3. map population density as a **choropleth**; 4. attach a country to every city with a **spatial join**, and fix the cities that miss; 5. **dissolve** countries into continents and save everything to a GeoPackage. Every number and figure below came from running [`geopandas_basics.py`](https://www.technicalgis.com/code/geopandas-getting-started/). ## Set up the environment Create a fresh virtual environment and install the pinned versions we tested with: `Terminal` ```bash python -m venv .venv .venv\Scripts\activate # macOS/Linux: source .venv/bin/activate pip install geopandas==0.13.2 shapely==2.0.7 pyogrio==0.9.0 pandas==2.0.3 matplotlib==3.7.5 mapclassify==2.6.0 ``` `mapclassify` gives `GeoDataFrame.plot()` its classification schemes (quantiles, natural breaks and so on). `pyogrio` is the fast GDAL-based reader that GeoPandas 1.x now uses by default. ## Step 1: read the data GeoPandas can read a zipped shapefile directly, so there is no need to unzip it first. The script downloads two Natural Earth 1:50m layers once and caches them in `./data`: `geopandas_basics.py` ```python import geopandas as gpd countries = gpd.read_file(fetch("ne_50m_admin_0_countries.zip")) places = gpd.read_file(fetch("ne_50m_populated_places.zip")) print(type(countries)) print(countries.shape) print(countries.crs) print(countries.geom_type.value_counts()) ``` `Output` ```text (242, 169) EPSG:4326 Polygon 123 MultiPolygon 119 Name: count, dtype: int64 ``` There are 242 countries and 169 attribute columns. The CRS is `EPSG:4326`, plain longitude and latitude in degrees. Almost half the features are `MultiPolygon`s: any country with islands or exclaves is stored as several parts. Most of those 169 columns are labels in other languages, so keep only what we need and give the columns friendlier names: ```python countries = countries[["NAME", "ISO_A3_EH", "CONTINENT", "POP_EST", "GDP_MD", "geometry"]] countries = countries.rename(columns={"NAME": "name", "ISO_A3_EH": "iso3", "CONTINENT": "continent", "POP_EST": "pop_est", "GDP_MD": "gdp_musd"}) print(countries.head()) ``` `Output` ```text name iso3 continent pop_est gdp_musd geometry 0 Zimbabwe ZWE Africa 14645468.0 21440 POLYGON ((31.28789 -22.40205, 31.19727 -22.344... 1 Zambia ZMB Africa 17861030.0 23309 POLYGON ((30.39609 -15.64307, 30.25068 -15.643... 2 Yemen YEM Asia 29161922.0 22581 MULTIPOLYGON (((53.08564 16.64839, 52.58145 16... 3 Vietnam VNM Asia 96462106.0 261921 MULTIPOLYGON (((104.06396 10.39082, 104.08301 ... 4 Venezuela VEN South America 28515829.0 482359 MULTIPOLYGON (((-60.82119 9.13838, -60.94141 9... ``` ## Step 2: reproject before you measure `countries.area` in EPSG:4326 would return square _degrees_, which mean nothing, because a degree of longitude shrinks towards the poles. To measure areas, reproject to an **equal-area** projection first. World Mollweide (`ESRI:54009`) works well for global data: ```python countries_ea = countries.to_crs("ESRI:54009") # World Mollweide (equal area) countries_ea["area_km2"] = countries_ea.area / 1e6 countries_ea["density"] = countries_ea["pop_est"] / countries_ea["area_km2"] print(countries_ea.nlargest(5, "area_km2")[["name", "area_km2", "pop_est", "density"]].round(1).to_string(index=False)) ``` `Output` ```text name area_km2 pop_est density Russia 16913027.1 144373535.0 8.5 Antarctica 12126825.6 4490.0 0.0 Canada 9894148.9 37589262.0 3.8 United States of America 9460566.8 328239523.0 34.7 China 9392194.5 1397715000.0 148.8 ``` These areas are within a few percent of the official figures; Russia’s official area is about 17.1 million km². The 1:50m coastlines are generalised, so small islands and intricate coasts lose a little area. ![Two world maps side by side: on the left in geographic latitude/longitude, where Greenland and Antarctica look huge; on the right in Mollweide, where land areas are shown in true proportion.](https://www.technicalgis.com/_astro/01-projections.CBfTQdfo_ZCU9Aq.webp) _The same GeoDataFrame before and after to_crs(). Plate carrée (left) inflates high latitudes; Mollweide (right) preserves area, so area-based maths is valid._ ## Step 3: a choropleth in one call With `mapclassify` installed, `plot()` can classify a column and draw a legend in one call. Quantile classes suit skewed data like population density, which runs from almost zero in Greenland to tens of thousands per km² in city-states: ```python world = countries_ea[countries_ea["continent"] != "Antarctica"] fig, ax = plt.subplots(figsize=(10, 5)) world.plot(column="density", scheme="Quantiles", k=6, cmap="YlGnBu", legend=True, edgecolor="white", linewidth=0.25, ax=ax, legend_kwds={"title": "People per km²", "loc": "lower left", "fmt": "{:.0f}", "fontsize": 7}) ax.set_title("Population density by country (quantile classes)") ax.set_axis_off() fig.savefig("output/02-density-choropleth.png") ``` ![Choropleth world map in Mollweide projection. Dense countries such as India, Bangladesh, the UK and much of western Europe are dark blue; Canada, Australia, Russia and Mongolia are pale yellow.](https://www.technicalgis.com/_astro/02-density-choropleth.DsGDa3UN_1zkAtY.webp) _Population density in six quantile classes. Each class holds roughly the same number of countries, so the colours spread evenly despite the skewed data._ ## Step 4: spatial joins — which country is each city in? The populated places layer has 1,251 points but no reliable country column. `gpd.sjoin()` adds one by testing every point against every polygon. GeoPandas uses a spatial index, so this takes milliseconds rather than the 300,000 brute-force tests it would otherwise need. Both layers must use the same CRS: ```python places = places[["NAME", "POP_MAX", "geometry"]].rename(columns={"NAME": "city", "POP_MAX": "city_pop"}) places_ea = places.to_crs(countries_ea.crs) joined = gpd.sjoin(places_ea, countries_ea[["name", "continent", "geometry"]], how="inner", predicate="within") print(joined[["city", "city_pop", "name", "continent"]].head().to_string()) print(len(places_ea) - len(joined), "of", len(places_ea), "places fall outside every country polygon") ``` `Output` ```text city city_pop name continent 0 Bombo 75000 Uganda Africa 1 Fort Portal 42670 Uganda Africa 13 Jinja 301619 Uganda Africa 14 Arua 250000 Uganda Africa 15 Mbale 402368 Uganda Africa 94 of 1251 places fall outside every country polygon ``` In a notebook the joined GeoDataFrame renders as a table. Columns from the right-hand layer (`name`, `continent`) now sit next to each city, and the index is preserved, which is why it jumps from 1 to 13: ![Notebook-style rendering of the first eight rows of the joined GeoDataFrame, showing city, city_pop, name, continent and projected POINT geometry columns for cities in Uganda.](https://www.technicalgis.com/_astro/05-joined-table.TOIjI634_Z6ABXY.webp) _joined.head(8) rendered as HTML, exactly as Jupyter would show it. Geometries are now in Mollweide metres._ ### Fixing the points that miss Seven percent of places matched nothing. Look at which ones and the pattern is clear: they are **coastal cities**. Natural Earth’s 1:50m coastline is simplified, so a city on the shore can sit a few hundred metres out to sea. `sjoin_nearest()` handles this. It assigns each leftover point to the closest polygon within a distance limit and reports how far away it was: ```python unmatched = places_ea[~places_ea.index.isin(joined.index)] nearest = gpd.sjoin_nearest(unmatched, countries_ea[["name", "continent", "geometry"]], max_distance=25_000, distance_col="dist_m") print(len(nearest), "recovered within 25 km") print(nearest.nlargest(5, "city_pop")[["city", "city_pop", "name", "dist_m"]].round(0).to_string()) joined = pd.concat([joined, nearest.drop(columns="dist_m")]) ``` `Output` ```text 92 recovered within 25 km city city_pop name dist_m 1226 New York 19040000 United States of America 1478.0 1228 Istanbul 10061000 Turkey 869.0 1233 Lagos 9466000 Nigeria 1576.0 1250 Hong Kong 7206000 Hong Kong 1327.0 1171 Salvador 3484000 Brazil 927.0 ``` New York, Istanbul, Lagos and Hong Kong, some of the largest cities on Earth, would have silently vanished from a simple `within` join, yet each lies less than 2 km from the generalised coastline. **Always count your unmatched rows after a spatial join.** Now we can aggregate as in plain pandas, for example to count the cities with more than a million people in each country: ```python big = joined[joined["city_pop"] >= 1_000_000] per_country = big.groupby("name").agg(cities_1m=("city", "count"), city_pop=("city_pop", "sum")) print(per_country.sort_values("cities_1m", ascending=False).head(8)) ``` `Output` ```text cities_1m city_pop name China 67 198833000 India 42 140258000 United States of America 40 132086838 Brazil 18 71900406 Russia 11 25533000 Japan 8 61636000 Nigeria 7 20462000 Indonesia 7 20886000 ``` Use the `.cx[]` coordinate indexer to clip to a bounding box and map Europe’s cities, with marker size scaled by population: ![Map of Europe in Mollweide projection with teal circles for cities sized by population; London, Paris, Madrid, Moscow, Istanbul and others are labelled.](https://www.technicalgis.com/_astro/03-spatial-join-europe.D5Mirk36_dpwni.webp) _European cities after the spatial join, sized by population. The ten largest are labelled._ ## Step 5: dissolve countries into continents `dissolve()` is GeoPandas’ _group by_ for geometry. It merges polygons that share a value and aggregates the other columns: ```python continents = countries_ea[countries_ea["continent"] != "Seven seas (open ocean)"].dissolve( by="continent", aggfunc={"pop_est": "sum", "gdp_musd": "sum", "area_km2": "sum"}) continents["density"] = continents["pop_est"] / continents["area_km2"] print(continents.drop(columns="geometry").sort_values("pop_est", ascending=False)) ``` `Output` ```text pop_est gdp_musd area_km2 density continent Asia 4,565,778,126.0 33555734 31,231,657.3 146.2 Africa 1,307,986,092.3 2459078 30,053,513.8 43.5 Europe 746,437,653.0 21642328 22,900,334.1 32.6 North America 585,351,846.0 25115795 24,230,324.1 24.2 South America 427,066,661.0 3852015 17,702,395.3 24.1 Oceania 42,430,881.0 1663441 8,539,629.5 5.0 Antarctica 4,490.0 898 12,126,825.6 0.0 ``` Europe’s 22.9 million km² is much larger than you might expect because Natural Earth assigns all of Russia, including Siberia, to Europe. Check how a dataset defines its groups before you quote aggregated figures. ![World map in Mollweide projection with each continent filled in a single pastel colour and internal country borders removed.](https://www.technicalgis.com/_astro/04-dissolve-continents.DeOChQPQ_Z27Q7oJ.webp) _After dissolve(): one feature per continent, with internal borders removed._ ## Step 6: save the results GeoPackage is a single-file, open-standard format that holds several layers, with none of the shapefile’s 10-character column-name limit: ```python gpkg = OUT / "natural_earth_analysis.gpkg" countries_ea.to_file(gpkg, layer="countries", driver="GPKG") continents.to_file(gpkg, layer="continents", driver="GPKG") joined.drop(columns="index_right").to_file(gpkg, layer="places_joined", driver="GPKG") ``` `Output` ```text natural_earth_analysis.gpkg 3.12 MB ``` Drag the `.gpkg` into QGIS and you will see all three layers, already in Mollweide. ## What you learned - A `GeoDataFrame` is a pandas DataFrame plus a geometry column and a CRS. - **Reproject before measuring**: use an equal-area CRS for areas and a local projected CRS for distances. - `sjoin()` joins by location; count the unmatched rows and use `sjoin_nearest()` for points that fall just outside polygons. - `dissolve()` is a spatial `groupby`. - `plot(column=..., scheme=...)` gives you a classified choropleth in one line. ## Next steps - Speed up joins and overlays on millions of features in [Fast vector geoprocessing with Shapely 2 and GeoParquet](https://www.technicalgis.com/tutorials/fast-vector-geoprocessing/). - Turn these GeoDataFrames into interactive maps in [Interactive web maps with Folium](https://www.technicalgis.com/tutorials/interactive-web-maps-folium/). - Find national boundary data for almost every country in our [GIS data sources catalogue](https://www.technicalgis.com/data-sources/). ## Data and licences - Admin 0 countries and populated places, 1:50m: [Natural Earth](https://www.naturalearthdata.com/), public domain. Downloaded from `naciscdn.org` in October 2026; population estimates are those bundled with the dataset. - Code: [`geopandas_basics.py`](https://www.technicalgis.com/code/geopandas-getting-started/), tested with Python 3.8.10, GeoPandas 0.13.2, Shapely 2.0.7 and pandas 2.0.3. Newer GeoPandas 1.x releases run it unchanged. ## Key terms - **CRS**: A coordinate reference system defines how coordinates map to places on Earth, for example EPSG:4326 (longitude/latitude) or ESRI:54009 (World Mollweide). - **Equal-area projection**: A map projection that preserves the relative size of areas, so area calculations are meaningful. - **Spatial join**: A join that matches rows by the spatial relationship of their geometries (within, intersects, nearest) rather than by a shared key. - **Choropleth**: A map that shades areas according to a classified attribute value, such as population density. ## Data for this tutorial - [Natural Earth](https://www.technicalgis.com/data-sources/#src-natural-earth) · Public-domain countries and populated places - [Boundary data for every country](https://www.technicalgis.com/data-sources/) Licences and formats for each source are listed in the [GIS data sources catalogue](https://www.technicalgis.com/data-sources/). --- Source: https://www.technicalgis.com/tutorials/geopandas-getting-started/