Plotting a GeoTIFF with Correct Extent and Axes

Pass the dataset bounds to imshow as an extent, in left-right-bottom-top order:

import matplotlib.pyplot as plt
import rasterio

with rasterio.open("scene.tif") as src:
    arr = src.read(1, masked=True)
    b = src.bounds                       # BoundingBox(left, bottom, right, top)

fig, ax = plt.subplots(figsize=(8, 8))
ax.imshow(arr, extent=(b.left, b.right, b.bottom, b.top), interpolation="nearest")
ax.ticklabel_format(style="plain")

The reordering on that third line is the whole trick, and it is where almost every flipped-raster bug comes from. This page belongs to rendering rasters with matplotlib in Visualization, Tiling & Web Delivery.

Why the Default Axes Are Wrong

imshow was designed for images, where row zero is the top and the axis unit is a pixel. A north-up raster also has its highest northing in row zero, so the image orientation happens to be right — but the axis labels are pixel indices, which tell a reader nothing about where they are looking.

Pixel axes against map axes With no extent the axes run from zero to the pixel dimensions and the vertical axis increases downward, so a reader cannot relate the figure to a map. With the bounds passed as an extent the axes read as eastings and northings, increasing in the conventional directions, and any vector layer in the same CRS can be drawn on top. no extent extent from src.bounds 0 10980 0 10980 axis units are pixels; y increases downward 699960 809760 9900000 9790200 axis units are metres; vectors overlay directly Same array, same orientation — only the axes changed, and with them everything you can do next.

The practical payoff is not the labels themselves but what they enable. Once the axes are in map units, a vector layer plots on top with no transformation, a scale bar can be drawn in data coordinates, and a reader can compare two figures of adjacent scenes.

Environment & Setup

Package Version pin Used for
rasterio >=1.3.0 Dataset bounds, CRS and masked reads
matplotlib >=3.8 imshow, tick formatting, aspect control
geopandas >=0.14 Optional vector overlay 
pip install "rasterio>=1.3.0" "matplotlib>=3.8" "geopandas>=0.14"

Complete Working Example

import matplotlib.pyplot as plt
import numpy as np
import rasterio


def plot_geotiff(path: str, band: int = 1, *, target_px: int = 1600,
                 cmap: str = "viridis", km_axes: bool = False):
    """Plot a band in map coordinates with sane axes."""
    with rasterio.open(path) as src:
        factor = max(1, int(max(src.width, src.height) / target_px))
        arr = src.read(band,
                       out_shape=(max(1, src.height // factor),
                                  max(1, src.width // factor)),
                       masked=True)
        b = src.bounds
        crs = src.crs

    # NOTE the reordering: BoundingBox is (left, bottom, right, top),
    # imshow wants (left, right, bottom, top)
    extent = (b.left, b.right, b.bottom, b.top)
    unit = crs.linear_units if crs and crs.is_projected else "degrees"
    if km_axes and crs and crs.is_projected:
        extent = tuple(v / 1000 for v in extent)
        unit = "km"

    valid = arr.compressed()
    vmin, vmax = np.percentile(valid, [2, 98]) if valid.size else (None, None)

    fig, ax = plt.subplots(figsize=(8, 8))
    im = ax.imshow(arr, extent=extent, cmap=cmap, vmin=vmin, vmax=vmax,
                   interpolation="nearest", origin="upper")
    ax.set_aspect("equal")                    # square pixels stay square
    ax.ticklabel_format(style="plain")        # no 1e6 offset
    ax.set_xlabel(f"Easting ({unit})")
    ax.set_ylabel(f"Northing ({unit})")
    fig.colorbar(im, ax=ax, shrink=0.75, pad=0.02)
    return fig, ax


if __name__ == "__main__":
    fig, ax = plot_geotiff("S2A_36NYF_20260614_B08.tif", km_axes=True)
    fig.savefig("b08.png", dpi=180, bbox_inches="tight")

origin="upper" is the default and is correct here: the array’s first row is the northernmost, and the extent’s top value corresponds to it. Setting origin="lower" without also swapping the extent produces an image that is upside down and mislabelled, which is worse than either alone.

Overlaying Vector Data

One axes object, two data sources, one CRS The raster defines the axes in its own projected CRS. A vector layer arriving in WGS84 must be reprojected into that CRS before plotting; once it is, both draw onto the same axes with no further transformation and the geometries land exactly where they belong. The axes belong to the raster; everything else must join it raster, EPSG:32636 defines the extent and the axis units vector, EPSG:4326 degrees — cannot be drawn on these axes to_crs(src.crs) now in metres, on the same grid as the image one axes object gdf.plot(ax=ax, ...) lands exactly right 
import geopandas as gpd

fields = gpd.read_file("field_boundaries.gpkg")

with rasterio.open("scene.tif") as src:
    fields = fields.to_crs(src.crs)          # essential, and easy to forget

fields.boundary.plot(ax=ax, edgecolor="#b71c1c", linewidth=0.8)
ax.set_xlim(extent[0], extent[1])            # keep the raster's extent
ax.set_ylim(extent[2], extent[3])

Resetting the limits after plotting the vector layer matters when the vector data extends beyond the scene: geopandas expands the axes to fit its own data, which silently zooms the figure out and leaves the raster as a small square in the corner.

Verification

import rasterio

with rasterio.open("scene.tif") as src:
    b = src.bounds
    assert b.left < b.right and b.bottom < b.top, "unexpected bounds ordering"
    assert src.transform.e < 0, "not a north-up raster; extent logic differs"

extent = (b.left, b.right, b.bottom, b.top)
assert extent[0] < extent[1] and extent[2] < extent[3]
The two transform elements that decide orientation The affine transform's a element is the pixel width and its e element is the pixel height, which is negative for a north-up raster because row index increases southward. A positive e means a south-up file, where the same extent and origin settings produce an inverted image. transform.e is the sign that matters north-up (almost always) a = +10, e = -10 row 0 is the northernmost origin="upper" is correct south-up (rare) a = +10, e = +10 row 0 is the southernmost needs origin="lower" Assert the sign rather than assuming it; the failure is a silently mirrored map.

The transform.e < 0 assertion is the one worth keeping. A north-up raster has a negative y pixel size, and the rare south-up file needs the opposite origin — silently plotting one upside down is the kind of error that survives a whole report. A quick visual confirmation is to overlay a coastline or a road layer: if it matches, the geometry is right; if it is mirrored about the horizontal axis, the extent’s middle values were swapped. Rotated transforms — where the b and d elements are non-zero — break the extent approach entirely, because imshow can only draw an axis-aligned rectangle; those files must be reprojected onto a north-up grid before they can be plotted this way.

Common Errors

The image is flipped vertically

src.bounds was passed straight through as the extent. Reorder it to (left, right, bottom, top). Unpacking the bounding box into named variables rather than positionally makes the mistake much harder to repeat.

The figure is stretched into a rectangle

The aspect ratio defaults to filling the axes. Call ax.set_aspect("equal") so square ground pixels stay square. If the pixels are genuinely non-square — some resampled products are — pass the ratio of the transform’s a and e magnitudes instead of "equal".

The vector overlay zooms the figure out

geopandas expanded the axes to fit geometry extending beyond the scene. Reset set_xlim and set_ylim to the raster extent after plotting the vector layer, or clip the layer to the raster bounds first.

Axis labels read +3.99e5

Matplotlib’s offset notation. Disable it with ax.ticklabel_format(style="plain"), or work in kilometres.

Frequently Asked Questions

Q: What order does the extent tuple use? Left, right, bottom, top — which is not the order rasterio’s bounds returns them in. BoundingBox gives left, bottom, right and top, so the two middle values must be swapped. Getting it wrong flips the image vertically without any error.

Q: Why do the axes show a scientific offset? Matplotlib factors out a common offset when coordinate values are large, which UTM eastings and northings always are. Call ax.ticklabel_format(style=‘plain’) to disable it, or divide the extent by a thousand and label the axes in kilometres.

Q: How do I overlay a shapefile on the raster? Reproject the vector layer to the raster’s CRS and plot it onto the same axes. Because the axes are already in map coordinates, the geometries land in the right place with no further transformation.

Q: Does this work for a raster in geographic coordinates? Yes, with one caveat: an equal aspect ratio in degrees stretches the map away from the equator, because a degree of longitude is shorter than a degree of latitude there. Either accept the distortion, set the aspect to the secant of the mean latitude, or reproject to a projected CRS before plotting.