Plotting choropleth rasters
import rioxarray as rxr
import matplotlib.pyplot as plt
from cartopy.crs import PlateCarree
import geomappy as mp
import geomappy.plot_utils
import os
os.chdir("../../../")
r = rxr.open_rasterio("data/wtd.tif", mask_and_scale=True)
a = r.values[0]
A contains a 2D raster of water table depth (Fan et al., 2017). To visualise this one can simply use matplotlib directly.
plt.imshow(a, cmap="Blues_r", vmax=100)
plt.colorbar()
plt.show()
At its simplests, geomappy does exactly this (although some esthetic
differences can be seen):
mp.plot_raster(a, cmap="Blues_r", vmax=100)
plt.show()
The biggest difference in workflow comes with the binning of the data. This gets handled internally instead of having to specify this with numpy outisde the plotting interface.
mp.plot_raster(a, bins=[0,0.1,0.5,1,2,5,10,25], cmap="Blues_r")
plt.show()
In this case, the colorbar can be converted into a true legend. For visibility we expand the figure size, by creating a new figure first.
f, ax = plt.subplots(figsize=(12, 12))
mp.plot_raster(a, bins=[0,0.1,0.5,1,2,5,10,25], cmap="Blues_r", legend='legend', ax=ax)
plt.show()
With a basemap
The functionality described above and in the section on basemaps can be applied here. First the bounds need to be extracted from the raster.
bounds = r.rio.bounds()
extent = bounds[0], bounds[2], bounds[1], bounds[3]
extent
(109.999999342, 155.000000419, -44.999998545, -8.999999499)
Then a basemap needs to be created
f, ax = plt.subplots(figsize=(6, 6), subplot_kw={'projection': PlateCarree()})
ax.coastlines()
ax.set_extent(extent)
geomappy.plot_utils.add_ticks(ax, 10)
geomappy.plot_utils.add_gridlines(ax, 10)
<cartopy.mpl.gridliner.Gridliner at 0x18305f980>
Then this GeoAxes object needs to be passed to the plotting function.
f, ax = plt.subplots(figsize=(6, 6), subplot_kw={'projection': PlateCarree()})
ax.coastlines()
ax.set_extent(extent)
geomappy.plot_utils.add_ticks(ax, 10)
geomappy.plot_utils.add_gridlines(ax, 10)
mp.plot_raster(a, ax=ax, cmap="Blues_r", bins=[0,0.1,0.5,1,2,5,10,25], extent=extent)
plt.show()
