Python:使用多边形在给定的二维网格上创建蒙版

Question

我有一些多边形(加拿大各省)，用GeoPandas读入，想用它们创建一个掩码以应用于二维网格数据纬度-经度网格(使用 iris 从 netcdf 文件读取)。最终目标是只保留给定省份的数据，其余数据被屏蔽掉。因此掩码对于省内的网格框为 1，对于省外的网格框为 0 或 NaN。

---

可以从此处的 shapefile 中获取多边形: https://www.dropbox.com/s/o5elu01fetwnobx/CAN_adm1.shp?dl=0

我使用的netcdf文件可以在这里下载: https://www.dropbox.com/s/kxb2v2rq17m7lp7/t2m.20090815.nc?dl=0

---

我想这里有两种方法，但我正在为这两种方法而苦苦挣扎:

1) 使用多边形在经纬度网格上创建一个掩码，以便它可以应用于 python 之外的大量数据文件(首选)

2) 使用多边形遮蔽读入的数据，只提取感兴趣省份内的数据，进行交互处理。

到目前为止我的代码:

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!

非常感谢任何帮助或建议。

---

更新:按照下面@DPeterK 的出色解决方案，我的原始数据可以被屏蔽，给出以下内容:

Answer 1

看起来你已经开始了!从 shapefile 加载的几何图形公开了各种地理空间比较方法，在这种情况下，您需要 contains 方法。您可以使用它来测试立方体水平网格中的每个点是否包含在不列颠哥伦比亚省几何中。 (请注意，这不是一个快速操作!)您可以使用此比较来构建一个 2D 掩码数组，它可以应用于您的立方体的数据或以其他方式使用。

我已经编写了一个 Python 函数来执行上述操作 – 它需要一个立方体和一个几何图形，并为立方体的(指定)水平坐标生成一个 mask ，并将 mask 应用于立方体的数据。函数如下:

def geom_to_masked_cube(cube, geometry, x_coord, y_coord,
                        mask_excludes=False):
    """
    Convert a shapefile geometry into a mask for a cube's data.

    Args:

    * cube:
        The cube to mask.
    * geometry:
        A geometry from a shapefile to define a mask.
    * x_coord: (str or coord)
        A reference to a coord describing the cube's x-axis.
    * y_coord: (str or coord)
        A reference to a coord describing the cube's y-axis.

    Kwargs:

    * mask_excludes: (bool, default False)
        If False, the mask will exclude the area of the geometry from the
        cube's data. If True, the mask will include *only* the area of the
        geometry in the cube's data.

    .. note::
        This function does *not* preserve lazy cube data.

    """
    # Get horizontal coords for masking purposes.
    lats = cube.coord(y_coord).points
    lons = cube.coord(x_coord).points
    lon2d, lat2d = np.meshgrid(lons,lats)

    # Reshape to 1D for easier iteration.
    lon2 = lon2d.reshape(-1)
    lat2 = lat2d.reshape(-1)

    mask = []
    # Iterate through all horizontal points in cube, and
    # check for containment within the specified geometry.
    for lat, lon in zip(lat2, lon2):
        this_point = gpd.geoseries.Point(lon, lat)
        res = geometry.contains(this_point)
        mask.append(res.values[0])

    mask = np.array(mask).reshape(lon2d.shape)
    if mask_excludes:
        # Invert the mask if we want to include the geometry's area.
        mask = ~mask
    # Make sure the mask is the same shape as the cube.
    dim_map = (cube.coord_dims(y_coord)[0],
               cube.coord_dims(x_coord)[0])
    cube_mask = iris.util.broadcast_to_shape(mask, cube.shape, dim_map)

    # Apply the mask to the cube's data.
    data = cube.data
    masked_data = np.ma.masked_array(data, cube_mask)
    cube.data = masked_data
    return cube

如果您只需要 2D 蒙版，您可以在上述函数将其应用于立方体之前返回它。

要在您的原始代码中使用此函数，请在代码末尾添加以下内容:

geometry = BritishColumbia.geometry
masked_cube = geom_to_masked_cube(cube, geometry,
                                  'longitude', 'latitude',
                                  mask_excludes=True)

如果这没有掩盖任何东西，则很可能意味着您的立方体和几何体是在不同的范围内定义的。也就是说，立方体的经度坐标范围为 0°–360°，如果几何体的经度值范围为 -180°–180°，则包含测试将永远不会返回 True。您可以通过以下方式更改立方体的范围来解决此问题:

cube = cube.intersection(longitude=(-180, 180))