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].

enter image description here

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 shape files 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 inital 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?


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


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.


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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.