Extract data from NetCDF using a shapefile polygon in Python [duplicate]

Question

I need to subset data from a NetCDF using a specific shapefile. The data are sea surface temperature and ocean color at 1/4 degree resolution. I have 4 polygons describing the US. Northeast continental shelf large marine ecosystem and it's sub-components that I need to use to extract the data. I am working with monthly composite files from 1982-2014, so this data extraction routine needs to be automated. The files are already subsetted to the approximate working area grid of [35, 45, -80, -60].

Previously, we were converting HDF5 data files to rasters in R and processing them this way, but this method is really inefficient and I am sure there is a better solution in Python using the current NetCDF files.

Thus far I have been using GDAL and Fiona to read in the shapefiles and NetCDF4 to load the data files. I am not sure how to go about subsetting the data. I found this:

GDAL for Python: extracting subdomains from NetCDF file?

But I don't have the foggiest idea about how to subset a NetCDF file using anything other than a simple bounding box, which these polygons most certainly are not.

Point in polygon routines would probably take an eternity to work, but maybe I could subset the data using a smaller bounding box that is rotated to fit these shapes like this as an initial starting point and then do a point-in-poly search:

Subsetting a curvilinear netCDF file (ROMS model output) using a lon/lat bounding box.

Any ideas?

EDIT 1:

I just came across the OpenClimateGIS package which seems that it may fit the bill perfectly... I will have a go with this to see if I can get it to work: http://ncpp.github.io/ocgis/examples.html#advanced-subsetting

Answer 1

This might be adaptable for your needs.

If you don't mind calling the command line from python, you could do something like gdalwarp -cutline clip.shp -cl clip -crop_to_cutline input_raster output_raster_clipped.tif. -cwhere and -csql might be more appropriate gdalwarp options for selecting one of the four polygons for clipping.

Answer 2

Take a look at this: https://stackoverflow.com/questions/34585582/how-to-mask-the-specific-array-data-based-on-the-shapefile

What you want to keep in mind is that once you've loaded your NetCDF, you're working with a NumPy array.

What are you trying to output? Summary stats based on the polygon areas?

Anyhow, here's what I'd do:

1. Load your shapefile and get your areas into a supported format (aiming for the matplotlib mask process in the above link sounds good)
2. Load your NetCDF file and get the data into a single X, Y, T numpy array
3. Mask that array using the polygons (one at a time?)
4. Export your summary stats.

Answer 3

You can use rioxarray. Here is an example: https://corteva.github.io/rioxarray/stable/examples/clip_geom.html

import rioxarray
import geopandas

geodf = geopandas.read_file(...)
xds = rioxarray.open_rasterio(...)
clipped = xds.rio.clip(geodf.geometry.apply(mapping), geodf.crs)