Parsing rasterio’s geocoordinates

Converting a projection’s cartesian coordinates into 2D longitudes and latitudes.

These new coordinates might be handy for plotting and indexing, but it should be kept in mind that a grid which is regular in projection coordinates will likely be irregular in lon/lat. It is often recommended to work in the data’s original map projection (see imshow() and map projections).

import as ccrs
import matplotlib.pyplot as plt
import numpy as np
from rasterio.warp import transform

import xarray as xr

# Read the data
url = ""
da = xr.open_rasterio(url)

# Compute the lon/lat coordinates with rasterio.warp.transform
ny, nx = len(da["y"]), len(da["x"])
x, y = np.meshgrid(da["x"], da["y"])

# Rasterio works with 1D arrays
lon, lat = transform(, {"init": "EPSG:4326"}, x.flatten(), y.flatten())
lon = np.asarray(lon).reshape((ny, nx))
lat = np.asarray(lat).reshape((ny, nx))
da.coords["lon"] = (("y", "x"), lon)
da.coords["lat"] = (("y", "x"), lat)

# Compute a greyscale out of the rgb image
greyscale = da.mean(dim="band")

# Plot on a map
ax = plt.subplot(projection=ccrs.PlateCarree())
ax.coastlines("10m", color="r")

Total running time of the script: ( 0 minutes 1.971 seconds)

Gallery generated by Sphinx-Gallery