Geopandas plot shapefile on xarray with same legend

user_15185911 0
Viewed 308

I'm trying to create some maps of precipitation data (xarray) with a shapefile of the region of interest on top. However, when Python plots the figures, I get two seperate figures:

The precipitation data

The region of interest

When I open the data in QGIS they do appear on top of each other, so the coordinate systems do check out. Then I have an additional bonus question: I have to create multiple precipitation maps, on for a visual analysis it would be ideal if I could have the same legend (thus the same min/max for the colorbar) for each map. Anyone an idea how to proceed further?

My code so far:

def chirps_to_map(input1, input2, title):

    projection = input1 + input2

    plt.figure(figsize=(9, 9))

    projection['pr'].plot()

    watershed.plot()

    plt.title(title)

    plt.show()

    plt.close()

    projection.to_netcdf(str(path)+str(title)+".nc")

    return projection
python-xarray geopandas matplotlib python
1 Answers

This is a case where it's simpler to use the Matplotlib object-oriented API.

A nice general workflow might be

fig, ax = plt.subplot()
gdf.plot(ax=ax)    # Plot the vector data on the subplot
raster.plot(ax=ax) # Plot the raster data on the same subplot

Example

First, we get some sample raster+vector data

import xarray as xr
import geopandas as gpd
import matplotlib.pyplot as plt

da = xr.tutorial.load_dataset('ROMS_example').zeta.isel(ocean_time=0)
gdf = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
usa = gdf.loc[gdf['name'].eq('United States of America')]

Next, we plot both of the data on the same AxesSubplot

fig, ax = plt.subplots(figsize=(15, 10))
da.plot.pcolormesh(x='lon_rho', y='lat_rho', ax=ax)
usa.plot(ax=ax, edgecolor='red', color='none')

# Focus on the raster extent
ax.set_xlim(-95, -87)
ax.set_ylim(26, 32)

enter image description here

Bonus: hvPlot way

hvPlot provides a nice unified API for interactive plotting with pandas, xarray, and many other libraries, and might be of interest to people stumbling upon this answer.

Plotting both vector and raster data is rather easy, simply use the * operator.

import hvplot.pandas
import hvplot.xarray

usa.hvplot(geo=True) * da.hvplot.quadmesh(x='lon_rho', y='lat_rho', geo=True)

user_17856610
Related