How to Plot Maps in Python with GeoPandas and Matplotlib

On this page

Problem statement

GeoDataFrame.plot() is the fastest check that a layer loaded where you expected, and the first step towards a static map for a report. This page covers it on shapefiles and GeoJSON:

check imported GIS data visually
plot points, lines, and polygons
create a simple thematic map
save the result to an image file

If you need a fast way to visualize vector spatial data in Python, this is the usual approach. Every example runs on real data for Rutland — OpenStreetMap places, roads and land use, and the ONS ward boundaries — plus the UK's local authority districts for the choropleth (set up under Code examples).

Quick answer

Two layers — boundary polygons and points — drawn together in the same CRS.

To plot spatial data with GeoPandas, read the file into a GeoDataFrame and call .plot():

import geopandas as gpd
import matplotlib.pyplot as plt

gdf = gpd.read_file("data/places.geojson")
gdf.plot()

plt.show()

This works well for:

quick static maps
exploratory GIS analysis
simple output for reports or automation scripts

For more control, pass styling arguments to .plot() and use Matplotlib to adjust the figure.

What you get — from Load spatial data into a GeoDataFrame

      WD24CD  ...                                           geometry
0  E05012632  ...  POLYGON ((486266.5 310507.598, 486329.797 3104...
1  E05012633  ...  POLYGON ((483986.402 307395.898, 484006.1 3073...
2  E05012634  ...  POLYGON ((489766.198 318469.798, 489812.098 31...

[3 rows x 3 columns]
EPSG:27700

Load, plot, style, layer and save a map

Understanding these five objects removes most of the guesswork from styling a map.

Load spatial data into a GeoDataFrame

GeoPandas can read common vector GIS formats including:

shapefile
GeoJSON
GeoPackage

Use gpd.read_file() to load the data:

import geopandas as gpd

gdf = gpd.read_file("data/admin_boundaries.shp")
print(gdf.head(3))
print(gdf.crs)

      WD24CD  ...                                           geometry
0  E05012632  ...  POLYGON ((486266.5 310507.598, 486329.797 3104...
1  E05012633  ...  POLYGON ((483986.402 307395.898, 484006.1 3073...
2  E05012634  ...  POLYGON ((489766.198 318469.798, 489812.098 31...

[3 rows x 3 columns]
EPSG:27700

To plot correctly, the data must contain a geometry column. You should also confirm that the file loaded with the expected coordinate reference system.

You can do the same with GeoJSON:

roads = gpd.read_file("data/roads.geojson")
print(len(roads), roads.geom_type.unique().tolist(), roads.crs)

7292 ['LineString'] EPSG:4326

Plot the map with GeoPandas

The simplest way to create a map is with the GeoDataFrame .plot() method:

import matplotlib.pyplot as plt

gdf.plot()
plt.show()

Default behavior depends on geometry type:

points are drawn as markers
lines are drawn as line paths
polygons are drawn as filled shapes

This is the fastest way to verify that the data loaded and that the geometry is where you expect it to be.

Customize the map appearance with Matplotlib

For clearer output, create a Matplotlib figure and axis first, then pass the axis to GeoPandas:

fig, ax = plt.subplots(figsize=(8, 6))

gdf.plot(
    ax=ax,
    color="lightblue",
    edgecolor="black",
    linewidth=0.8
)

ax.set_title("Administrative Boundaries")
plt.show()

Useful options include:

figsize for map size
color for fill or marker color
edgecolor for polygon borders
linewidth for line or border thickness
markersize for point layers

Plot a column as a thematic map

To create a choropleth map, pass a numeric column to column=:

districts = gpd.read_file("data/districts.geojson")
fig, ax = plt.subplots(figsize=(10, 8))

districts.plot(
    column="population",
    ax=ax,
    cmap="OrRd",
    legend=True,
    edgecolor="black",
    linewidth=0.5
)

ax.set_title("Population by District")
plt.show()

This is useful when you want to visualize an attribute such as:

population
area
density
land value
any other numeric field

For categorical fields, you can also use column=, though the output is usually clearer when the number of classes is small.

Plot multiple layers on the same map

Roads over ward boundaries, places over roads: layers drawn on one ax line up only if they share a CRS.

import geopandas as gpd
import matplotlib.pyplot as plt

boundaries = gpd.read_file("data/admin_boundaries.shp")
schools = gpd.read_file("data/schools.geojson")

schools = schools.to_crs(boundaries.crs)

fig, ax = plt.subplots(figsize=(10, 8))

boundaries.plot(
    ax=ax,
    color="white",
    edgecolor="black",
    linewidth=0.8
)

schools.plot(
    ax=ax,
    color="red",
    markersize=20
)

ax.set_title("Schools Within Administrative Boundaries")
ax.set_axis_off()
plt.show()

Use this pattern when you need to check alignment between layers or produce a simple reference map. The to_crs line is the one that matters. The schools come from OpenStreetMap in longitude and latitude and the wards from the ONS in British National Grid; drawn without it, the two layers share an axes and nothing else:

fig, ax = plt.subplots()
boundaries.plot(ax=ax)
gpd.read_file("data/schools.geojson").plot(ax=ax, color="red")
x0, x1 = ax.get_xlim()
print(f"x axis spans {x0:,.1f} to {x1:,.1f}: the schools sit at about -0.7, the wards at about 490,000")
plt.close(fig)

x axis spans -25,315.6 to 531,610.4: the schools sit at about -0.7, the wards at about 490,000

Remove axes and improve map readability

Map axes are often unnecessary in a simple static output. You can hide them:

parks = gpd.read_file("data/landuse.shp").query("fclass == 'park'")
fig, ax = plt.subplots(figsize=(8, 8))

parks.plot(ax=ax, color="lightgreen", edgecolor="gray")
ax.set_title("Park Boundaries")
ax.set_axis_off()

plt.show()

For thematic maps, keep legend=True if the attribute values need explanation.

Save the map to an image file

To export the map for a report or automated workflow, save the figure with Matplotlib:

from pathlib import Path
Path("output").mkdir(exist_ok=True)

fig, ax = plt.subplots(figsize=(10, 8))

districts.plot(
    column="population",
    ax=ax,
    cmap="Blues",
    legend=True,
    edgecolor="black"
)

ax.set_title("Population Map")
ax.set_axis_off()

plt.savefig("output/population_map.png", dpi=300, bbox_inches="tight")
plt.close(fig)
print(Image.open("output/population_map.png").size, "pixels")

(1360, 1971) pixels

This is the standard pattern for repeatable GIS reporting scripts.

Rutland wards, roads and schools as maps

The data used on this page

Geofabrik's OpenStreetMap extract for Rutland (its dated 1 January 2026 snapshot) — the places, roads, land use and schools layers, written as the GeoJSON and shapefiles the examples read — the ONS December 2024 ward boundaries for Rutland as admin_boundaries.shp, and the UK's local authority districts with their mid-2025 population from the ONS estimates on NOMIS as districts.geojson.

import io
import zipfile
from pathlib import Path

import geopandas as gpd
import matplotlib
import pandas as pd
import requests
from PIL import Image

matplotlib.use("Agg")
import matplotlib.pyplot as plt

data = Path("data")
data.mkdir(exist_ok=True)
if not (data / "landuse.shp").exists():
    reply = requests.get("https://download.geofabrik.de/europe/united-kingdom/england/"
                         "rutland-260101-free.shp.zip", timeout=300)
    reply.raise_for_status()
    zipfile.ZipFile(io.BytesIO(reply.content)).extractall("osm")
    places = gpd.read_file("osm/gis_osm_places_free_1.shp")
    places[places["fclass"].isin(["town", "village"])].to_file(data / "places.geojson")
    gpd.read_file("osm/gis_osm_roads_free_1.shp").to_file(data / "roads.geojson")
    gpd.read_file("osm/gis_osm_landuse_a_free_1.shp").to_file(data / "landuse.shp")
    pois = pd.concat([gpd.read_file(f"osm/gis_osm_{layer}_free_1.shp") for layer in ("pois", "pois_a")])
    schools = pois[pois["fclass"] == "school"].copy()
    schools["geometry"] = schools.to_crs(27700).representative_point().to_crs(4326)
    schools.to_file(data / "schools.geojson")

ONS = "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
if not (data / "admin_boundaries.shp").exists():
    reply = requests.get(ONS + "Wards_December_2024_Boundaries_UK_BGC/FeatureServer/0/query", timeout=120,
                         params={"where": "LAD24NM='Rutland'", "outFields": "WD24CD,WD24NM",
                                 "outSR": 27700, "f": "geojson"})
    reply.raise_for_status()
    gpd.read_file(io.BytesIO(reply.content)).to_file(data / "admin_boundaries.shp")
if not (data / "districts.geojson").exists():
    reply = requests.get(ONS + "Local_Authority_Districts_December_2024_Boundaries_UK_BUC/FeatureServer/0/query",
                         timeout=120, params={"where": "1=1", "outFields": "LAD24CD,LAD24NM",
                                              "outSR": 4326, "f": "geojson"})
    reply.raise_for_status()
    population = pd.read_csv("https://www.nomisweb.co.uk/api/v01/dataset/NM_2002_1.data.csv?geography=TYPE424"
                             "&date=latest&gender=0&c_age=200&measures=20100"
                             "&select=geography_code,obs_value")
    (gpd.read_file(io.BytesIO(reply.content))
        .merge(population.rename(columns={"OBS_VALUE": "population"}), left_on="LAD24CD", right_on="GEOGRAPHY_CODE")
        .dropna(subset=["population"]).drop(columns="GEOGRAPHY_CODE").to_file(data / "districts.geojson"))

for name in ["places.geojson", "roads.geojson", "landuse.shp", "schools.geojson", "admin_boundaries.shp",
             "districts.geojson"]:
    layer = gpd.read_file(data / name)
    print(f"{name:21} {len(layer):>6,} {layer.geom_type.iloc[0]:<12} {layer.crs}")

places.geojson            54 Point        EPSG:4326
roads.geojson          7,292 LineString   EPSG:4326
landuse.shp            2,405 Polygon      EPSG:4326
schools.geojson           28 Point        EPSG:4326
admin_boundaries.shp      15 Polygon      EPSG:27700
districts.geojson        361 Polygon      EPSG:4326

Example 1: Plot a shapefile as a basic map

Use this when you want to confirm that a shapefile loaded correctly.

import geopandas as gpd
import matplotlib.pyplot as plt

gdf = gpd.read_file("data/landuse.shp")

gdf.plot()
plt.show()

Example 2: Style polygon boundaries and fill colors

This is a better choice for a simple reporting map.

import geopandas as gpd
import matplotlib.pyplot as plt

gdf = gpd.read_file("data/landuse.shp")

fig, ax = plt.subplots(figsize=(9, 7))
gdf.plot(
    ax=ax,
    color="beige",
    edgecolor="darkslategray",
    linewidth=0.6
)

ax.set_title("Land Use Areas")
ax.set_axis_off()
plt.show()

Example 3: Create a choropleth map from an attribute column

Use a numeric field to color polygons.

import geopandas as gpd
import matplotlib.pyplot as plt

gdf = gpd.read_file("data/districts.geojson").to_crs(27700)
gdf["pop_density"] = gdf["population"] / (gdf.area / 1e6)          # people per km²

fig, ax = plt.subplots(figsize=(10, 8))
gdf.plot(
    column="pop_density",
    ax=ax,
    cmap="YlGnBu",
    legend=True,
    edgecolor="black",
    linewidth=0.4
)

ax.set_title("Population Density")
ax.set_axis_off()
plt.show()
print(gdf.nlargest(3, "pop_density")[["LAD24NM", "pop_density"]].round(0).to_string(index=False))
palest = (gdf["pop_density"] < gdf["pop_density"].max() * 0.1).mean()
print(f"{palest:.0%} of districts fall in the palest tenth of the colour ramp")

      LAD24NM  pop_density
Tower Hamlets      17328.0
    Islington      14824.0
      Hackney      13424.0
70% of districts fall in the palest tenth of the colour ramp

Example 4: Plot point data on a clean map layout

Point datasets often need marker styling to be readable.

import geopandas as gpd
import matplotlib.pyplot as plt

gdf = gpd.read_file("data/schools.geojson")

fig, ax = plt.subplots(figsize=(8, 8))
gdf.plot(
    ax=ax,
    color="red",
    markersize=20
)

ax.set_title("School Locations")
ax.set_axis_off()
plt.show()

Example 5: Save the plotted map to a PNG file

Use this pattern in automated GIS workflows.

import geopandas as gpd
import matplotlib.pyplot as plt

gdf = gpd.read_file("data/admin_boundaries.shp")

fig, ax = plt.subplots(figsize=(10, 8))
gdf.plot(ax=ax, color="lightgray", edgecolor="black", linewidth=0.5)

ax.set_title("Boundary Map")
ax.set_axis_off()

plt.savefig("output/boundary_map.png", dpi=300, bbox_inches="tight")
plt.close(fig)

How GeoPandas hands drawing to Matplotlib

GeoPandas and Matplotlib work together:

GeoPandas reads and stores spatial data in a GeoDataFrame
Matplotlib renders the figure and handles styling

Six small changes, and the same data starts to look deliberate.

When you call gdf.plot(), GeoPandas converts the geometry into a Matplotlib plot. In many GIS workflows, this is enough for quick visualization.

This method is best for static maps. It is not intended for interactive web maps with pan, zoom, popups, or tiled basemaps.

Geometry type affects the output automatically:

Point geometries become markers
LineString geometries become lines
Polygon geometries become filled areas with optional borders

For many tasks, GeoPandas plotting is the fastest way to inspect data before moving on to reprojection, clipping, spatial joins, or export.

CRS, invalid geometry, large layers, static output

Coordinate reference systems can affect map output

If a layer appears in the wrong place or looks distorted, check its CRS:

print(gdf.crs)

EPSG:27700

If you are plotting multiple layers together, make sure they use the same CRS before plotting:

roads = roads.to_crs(boundaries.crs)

For single-layer plotting, reprojection is only necessary if your workflow requires a different CRS.

Missing or invalid geometry can affect plotting

Null or invalid geometries may cause blank output, incomplete plots, unexpected rendering, or errors.

Check for missing geometry:

print(gdf.geometry.isna().sum())

Check validity:

print(gdf.geometry.is_valid.value_counts())

True    15
Name: count, dtype: int64

If necessary, filter or repair problematic rows before plotting.

Large datasets may plot slowly

Very large shapefiles or GeoJSON files can be slow to render. Common fixes include:

selecting fewer columns
filtering rows
plotting a sample
simplifying geometry before plotting

This is common with parcel, road, or building datasets.

GeoPandas plots static maps only

If you need interactive maps, GeoPandas .plot() is not the right tool. This page covers static plotting for analysis, reporting, and script-based GIS workflows.

Internal links

FAQ

How do I plot a shapefile in Python with GeoPandas?

Read the shapefile with gpd.read_file() and call .plot():

import geopandas as gpd
import matplotlib.pyplot as plt

gdf = gpd.read_file("data/admin_boundaries.shp")
gdf.plot()
plt.show()

Can GeoPandas plot GeoJSON files directly?

Yes. GeoPandas can read GeoJSON with read_file() the same way as a shapefile:

gdf = gpd.read_file("data/places.geojson")
gdf.plot()

How do I color a map by an attribute column in GeoPandas?

Use the column argument and optionally add a legend:

districts.plot(column="population", cmap="OrRd", legend=True)

This is the standard way to create a choropleth map in GeoPandas.

How do I save a GeoPandas plot as a PNG image?

Create the plot, then save it with Matplotlib:

fig, ax = plt.subplots()
gdf.plot(ax=ax)
plt.savefig("map.png", dpi=300, bbox_inches="tight")
plt.close(fig)

This is useful for reports, dashboards, and batch GIS scripts.

How do I add a basemap or background tiles behind my GeoPandas plot?

GeoPandas .plot() does not provide basemaps itself; use the contextily library, which adds tiled background imagery to a Matplotlib axis. Reproject your layer to EPSG:3857 (Web Mercator) first, since most tile providers expect that CRS, then call contextily.add_basemap(ax).

Why does my choropleth use one color for everything or look washed out?

This usually happens when the column values are dominated by a few extremes, so most features fall into one color band. In Example 3's density map, inner London's 17,000 people per km² stretch the ramp so far that 70% of districts share its palest tenth. Pass a classification scheme such as scheme="quantiles" (which requires the mapclassify package) to spread features across classes more evenly.

How do I control the colorbar or legend size and position?

For a continuous colorbar, pass legend_kwds={"shrink": 0.5} to shrink it, or build a divider with mpl_toolkits.axes_grid1.make_axes_locatable for precise placement. For categorical legends, legend_kwds={"loc": "lower right"} moves the legend box.

How do I plot only the outlines of polygons with no fill?

Set the fill to none and keep an edge color, for example gdf.plot(ax=ax, facecolor="none", edgecolor="black"). This is useful for drawing boundary layers on top of a filled or basemap layer without hiding what is underneath.