Read and Plot Geotiff in python (xarray, cartopy)

When you’re working with geographic data you have a few choices

1. Work exclusively with GIS software
2. Work exclusively with python (or your favorite programming language)
3. Clean up your data using a pre-processor (python, geopandas, xarray) and then prettify using GIS software (QGIS, etc). This is the recommended choice.

Sometimes you can be someone like me who is very very stubborn and you end up doing either 1 or 2. It will be hard, but sometimes it’s worth it.

For this story we’ll talk about choice #2.

Say you have a pretty geotiff file and you want to integrate it into your plots. You can easily load your geotiff using xarray.open_rasterio.

da = xr.open_rasterio(’D:/Data/SRTM_Ph/srtm_luzon_clipped.tif’)
transform = Affine.from_gdal(*da.attrs[’transform’]) # this is important to retain the geographic attributes from the file

Now you want to plot it. Easiest and fastest way to do it is using matplotlib.pyplot.imshow()

fig = plt.figure(figsize=(16,8))
ax = fig.add_subplot(111)
ax.imshow(da.variable.data[0])
plt.show()

There’s one problem with this though: it treats your geotiff as an image and removes geographic data such as longitude and latitude. Notice the x and y axis in the image above showing pixel numbers rather than longs and lats.

Calling cartopy doesn’t help because the data is in a different space (ie the x and y values are different than the lons/lats). Forcing a projection simply flips your image without adding any geographic data.

crs=ccrs.PlateCarree()
fig = plt.figure(figsize=(16,8))
ax = fig.add_subplot(111, projection=crs)
ax.imshow(da.variable.data[0])
plt.show()

One solution is to create a meshgrid based on your geotiff information and then you can use this as inputs to matplotlib’s contourf or pcolormesh

# Create meshgrid from geotiff
nx, ny = da.sizes['x'], da.sizes['y']
x, y = np.meshgrid(np.arange(nx), np.arange(ny)) * transform
# Plot!
fig = plt.figure(figsize=(16,8))
ax = fig.add_subplot(111, projection=crs)
ax.coastlines(resolution='10m', alpha=0.1)
ax.contourf(x, y, da.variable.data[0], cmap='Greys')
ax.set_extent([lon_min, lon_max, lat_min, lat_max])
plt.show()

Tada!

To make things fancier, you can label your grids and axes

# Plot!
fig = plt.figure(figsize=(16,8))
ax = fig.add_subplot(111, projection=crs)
ax.coastlines(resolution='10m', alpha=0.1)
ax.contourf(x, y, da.variable.data[0], cmap='Greys')
ax.set_extent([lon_min, lon_max, lat_min, lat_max])
# Grid and Labels
gl = ax.gridlines(crs=crs, draw_labels=True, alpha=0.5)
gl.xlabels_top = None
gl.ylabels_right = None
xgrid = np.arange(lon_min-0.5, lon_max+0.5, 1.)
ygrid = np.arange(lat_min, lat_max+1, 1.)
gl.xlocator = mticker.FixedLocator(xgrid.tolist())
gl.ylocator = mticker.FixedLocator(ygrid.tolist())
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.xlabel_style = {'size': 14, 'color': 'black'}
gl.ylabel_style = {'size': 14, 'color': 'black'}
plt.show()

Tada

Because I’m lazy I will just share the whole notebook here