I have some polygons (Canadian provinces), read in with
GeoPandas, and want to use these to create a mask to apply to gridded data on a 2-d latitude-longitude grid (read from a netcdf file using
iris). An end goal would be to only have data for a given province remaining, with the rest of the data masked out. So the mask would be 1's for grid boxes within the province, and 0's or NaN's for grid boxes outside the province.
The polygons can be obtained from the shapefile here: https://www.dropbox.com/s/o5elu01fetwnobx/CAN_adm1.shp?dl=0
The netcdf file I am using can be downloaded here: https://www.dropbox.com/s/kxb2v2rq17m7lp7/t2m.20090815.nc?dl=0
I imagine there are two approaches here but I am struggling with both:
1) Use the polygon to create a mask on the latitude-longitude grid so that this can be applied to lots of datafiles outside of python (preferred)
2) Use the polygon to mask the data that have been read in and extract only the data inside the province of interest, to work with interactively.
My code so far:
import iris import geopandas as gpd #read the shapefile and extract the polygon for a single province #(province names stored as variable 'NAME_1') Canada=gpd.read_file('CAN_adm1.shp') BritishColumbia=Canada[Canada['NAME_1'] == 'British Columbia'] #get the latitude-longitude grid from netcdf file cubelist=iris.load('t2m.20090815.nc') cube=cubelist lats=cube.coord('latitude').points lons=cube.coord('longitude').points #create 2d grid from lats and lons (may not be necessary?) [lon2d,lat2d]=np.meshgrid(lons,lats)