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[0]
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)
 
#HELP!

Thanks very much for any help or advice.

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[0]
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)