Python：使用多边形在给定的二维网格上创建蒙版
Python: using polygons to create a mask on a given 2d grid
我有一些 多边形（加拿大各省），用 GeoPandas 读入，想用它们创建一个掩码以应用于 2 上的网格化数据-d 纬度-经度网格（使用 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 的出色解决方案，我的原始数据可以被屏蔽，给出以下内容：
看来你开局不错！从 shapefile 加载的几何图形公开了各种地理空间比较方法，在这种情况下，您需要 contains 方法。您可以使用它来测试立方体水平网格中的每个点是否包含在不列颠哥伦比亚省几何中。 （请注意，这 不是 快速操作！）您可以使用此比较来构建 2D 掩码数组，它可以应用于立方体的数据或以其他方式使用。
我编写了一个 Python 函数来执行上述操作 – 它需要一个立方体和一个几何体，并为立方体的（指定）水平坐标生成一个遮罩，并将遮罩应用于立方体的数据。函数如下：
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 蒙版，您可以 return 在上述函数将其应用于立方体之前。
要在您的原始代码中使用此函数，请在代码末尾添加以下内容：
geometry = BritishColumbia.geometry
masked_cube = geom_to_masked_cube(cube, geometry,
'longitude', 'latitude',
mask_excludes=True)
如果这没有掩盖任何东西，则很可能意味着您的立方体和几何图形是在不同的范围内定义的。也就是说，你的立方体的经度坐标 运行s 从 0°–360°，如果几何的经度值 运行 从 -180°–180°，那么包含测试永远不会 return True。您可以通过以下方式更改立方体的范围来解决此问题：
cube = cube.intersection(longitude=(-180, 180))
我找到了上面@DPeterK 发布的优秀解决方案的替代解决方案，它产生了相同的结果。它使用 matplotlib.path 来测试点是否包含在 外部坐标 中，这些坐标由从形状文件加载的 几何体 描述。 我发布这个是因为这个方法比@DPeterK 给出的方法快 10 倍（2:23 分钟 vs 25:56 分钟）。 我不确定什么更好：优雅的解决方案，或快速的蛮力解决方案。也许一个人可以两者兼得？！
此方法的一个复杂问题是某些几何图形是 MultiPolygons - 即形状由几个较小的多边形组成（在这种情况下，不列颠哥伦比亚省包括西海岸，不能用不列颠哥伦比亚省大陆的坐标来描述 Polygon)。 MultiPolygon 没有外部坐标，但单个多边形有，因此需要单独处理这些多边形。我发现最巧妙的解决方案是使用从 GitHub (https://gist.github.com/mhweber/cf36bb4e09df9deee5eb54dc6be74d26) 复制的函数，其中 'explodes' MultiPolygons 到一个单独的多边形列表中，然后可以单独处理。
下面列出了工作代码以及我的文档。抱歉，这不是最优雅的代码 - 我对 Python 比较陌生，我确信有很多不必要的 loops/neater 方法来做事！
import numpy as np
import iris
import geopandas as gpd
from shapely.geometry import Point
import matplotlib.path as mpltPath
from shapely.geometry.polygon import Polygon
from shapely.geometry.multipolygon import MultiPolygon
#-----
#FIRST, read in the target data and latitude-longitude grid from netcdf file
cubelist=iris.load('t2m.20090815.minus180_180.nc')
cube=cubelist[0]
lats=cube.coord('latitude').points
lons=cube.coord('longitude').points
#create 2d grid from lats and lons
[lon2d,lat2d]=np.meshgrid(lons,lats)
#create a list of coordinates of all points within grid
points=[]
for latit in range(0,241):
for lonit in range(0,480):
point=(lon2d[latit,lonit],lat2d[latit,lonit])
points.append(point)
#turn into np array for later
points=np.array(points)
#get the cube data - useful for later
fld=np.squeeze(cube.data)
#create a mask array of zeros, same shape as fld, to be modified by
#the code below
mask=np.zeros_like(fld)
#NOW, read the shapefile and extract the polygon for a single province
#(province names stored as variable 'NAME_1')
Canada=gpd.read_file('/Users/ianashpole/Computing/getting_province_outlines/CAN_adm_shp/CAN_adm1.shp')
BritishColumbia=Canada[Canada['NAME_1'] == 'British Columbia']
#BritishColumbia.geometry.type reveals this to be a 'MultiPolygon'
#i.e. several (in this case, thousands...) if individual polygons.
#I ultimately want to get the exterior coordinates of the BritishColumbia
#polygon, but a MultiPolygon is a list of polygons and therefore has no
#exterior coordinates. There are probably many ways to progress from here,
#but the method I have stumbled upon is to 'explode' the multipolygon into
#it's individual polygons and treat each individually. The function below
#to 'explode' the MultiPolygon was found here:
#https://gist.github.com/mhweber/cf36bb4e09df9deee5eb54dc6be74d26
#---define function to explode MultiPolygons
def explode_polygon(indata):
indf = indata
outdf = gpd.GeoDataFrame(columns=indf.columns)
for idx, row in indf.iterrows():
if type(row.geometry) == Polygon:
#note: now redundant, but function originally worked on
#a shapefile which could have combinations of individual polygons
#and MultiPolygons
outdf = outdf.append(row,ignore_index=True)
if type(row.geometry) == MultiPolygon:
multdf = gpd.GeoDataFrame(columns=indf.columns)
recs = len(row.geometry)
multdf = multdf.append([row]*recs,ignore_index=True)
for geom in range(recs):
multdf.loc[geom,'geometry'] = row.geometry[geom]
outdf = outdf.append(multdf,ignore_index=True)
return outdf
#-------
#Explode the BritishColumbia MultiPolygon into its constituents
EBritishColumbia=explode_polygon(BritishColumbia)
#Loop over each individual polygon and get external coordinates
for index,row in EBritishColumbia.iterrows():
print 'working on polygon', index
mypolygon=[]
for pt in list(row['geometry'].exterior.coords):
print index,', ',pt
mypolygon.append(pt)
#See if any of the original grid points read from the netcdf file earlier
#lie within the exterior coordinates of this polygon
#pth.contains_points returns a boolean array (true/false), in the
#shape of 'points'
path=mpltPath.Path(mypolygon)
inside=path.contains_points(points)
#find the results in the array that were inside the polygon ('True')
#and set them to missing. First, must reshape the result of the search
#('points') so that it matches the mask & original data
#reshape the result to the main grid array
inside=np.array(inside).reshape(lon2d.shape)
i=np.where(inside == True)
mask[i]=1
print 'fininshed checking for points inside all polygons'
#mask now contains 0's for points that are not within British Columbia, and
#1's for points that are. FINALLY, use this to mask the original data
#(stored as 'fld')
i=np.where(mask == 0)
fld[i]=np.nan
#Done.
我有一些 多边形（加拿大各省），用 GeoPandas 读入，想用它们创建一个掩码以应用于 2 上的网格化数据-d 纬度-经度网格（使用 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 的出色解决方案，我的原始数据可以被屏蔽，给出以下内容：
看来你开局不错！从 shapefile 加载的几何图形公开了各种地理空间比较方法，在这种情况下，您需要 contains 方法。您可以使用它来测试立方体水平网格中的每个点是否包含在不列颠哥伦比亚省几何中。 （请注意，这 不是 快速操作！）您可以使用此比较来构建 2D 掩码数组，它可以应用于立方体的数据或以其他方式使用。
我编写了一个 Python 函数来执行上述操作 – 它需要一个立方体和一个几何体，并为立方体的（指定）水平坐标生成一个遮罩，并将遮罩应用于立方体的数据。函数如下：
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 蒙版，您可以 return 在上述函数将其应用于立方体之前。
要在您的原始代码中使用此函数，请在代码末尾添加以下内容：
geometry = BritishColumbia.geometry
masked_cube = geom_to_masked_cube(cube, geometry,
'longitude', 'latitude',
mask_excludes=True)
如果这没有掩盖任何东西，则很可能意味着您的立方体和几何图形是在不同的范围内定义的。也就是说，你的立方体的经度坐标 运行s 从 0°–360°，如果几何的经度值 运行 从 -180°–180°，那么包含测试永远不会 return True。您可以通过以下方式更改立方体的范围来解决此问题：
cube = cube.intersection(longitude=(-180, 180))
我找到了上面@DPeterK 发布的优秀解决方案的替代解决方案，它产生了相同的结果。它使用 matplotlib.path 来测试点是否包含在 外部坐标 中，这些坐标由从形状文件加载的 几何体 描述。 我发布这个是因为这个方法比@DPeterK 给出的方法快 10 倍（2:23 分钟 vs 25:56 分钟）。 我不确定什么更好：优雅的解决方案，或快速的蛮力解决方案。也许一个人可以两者兼得？！
此方法的一个复杂问题是某些几何图形是 MultiPolygons - 即形状由几个较小的多边形组成（在这种情况下，不列颠哥伦比亚省包括西海岸，不能用不列颠哥伦比亚省大陆的坐标来描述 Polygon)。 MultiPolygon 没有外部坐标，但单个多边形有，因此需要单独处理这些多边形。我发现最巧妙的解决方案是使用从 GitHub (https://gist.github.com/mhweber/cf36bb4e09df9deee5eb54dc6be74d26) 复制的函数，其中 'explodes' MultiPolygons 到一个单独的多边形列表中，然后可以单独处理。
下面列出了工作代码以及我的文档。抱歉，这不是最优雅的代码 - 我对 Python 比较陌生，我确信有很多不必要的 loops/neater 方法来做事！
import numpy as np
import iris
import geopandas as gpd
from shapely.geometry import Point
import matplotlib.path as mpltPath
from shapely.geometry.polygon import Polygon
from shapely.geometry.multipolygon import MultiPolygon
#-----
#FIRST, read in the target data and latitude-longitude grid from netcdf file
cubelist=iris.load('t2m.20090815.minus180_180.nc')
cube=cubelist[0]
lats=cube.coord('latitude').points
lons=cube.coord('longitude').points
#create 2d grid from lats and lons
[lon2d,lat2d]=np.meshgrid(lons,lats)
#create a list of coordinates of all points within grid
points=[]
for latit in range(0,241):
for lonit in range(0,480):
point=(lon2d[latit,lonit],lat2d[latit,lonit])
points.append(point)
#turn into np array for later
points=np.array(points)
#get the cube data - useful for later
fld=np.squeeze(cube.data)
#create a mask array of zeros, same shape as fld, to be modified by
#the code below
mask=np.zeros_like(fld)
#NOW, read the shapefile and extract the polygon for a single province
#(province names stored as variable 'NAME_1')
Canada=gpd.read_file('/Users/ianashpole/Computing/getting_province_outlines/CAN_adm_shp/CAN_adm1.shp')
BritishColumbia=Canada[Canada['NAME_1'] == 'British Columbia']
#BritishColumbia.geometry.type reveals this to be a 'MultiPolygon'
#i.e. several (in this case, thousands...) if individual polygons.
#I ultimately want to get the exterior coordinates of the BritishColumbia
#polygon, but a MultiPolygon is a list of polygons and therefore has no
#exterior coordinates. There are probably many ways to progress from here,
#but the method I have stumbled upon is to 'explode' the multipolygon into
#it's individual polygons and treat each individually. The function below
#to 'explode' the MultiPolygon was found here:
#https://gist.github.com/mhweber/cf36bb4e09df9deee5eb54dc6be74d26
#---define function to explode MultiPolygons
def explode_polygon(indata):
indf = indata
outdf = gpd.GeoDataFrame(columns=indf.columns)
for idx, row in indf.iterrows():
if type(row.geometry) == Polygon:
#note: now redundant, but function originally worked on
#a shapefile which could have combinations of individual polygons
#and MultiPolygons
outdf = outdf.append(row,ignore_index=True)
if type(row.geometry) == MultiPolygon:
multdf = gpd.GeoDataFrame(columns=indf.columns)
recs = len(row.geometry)
multdf = multdf.append([row]*recs,ignore_index=True)
for geom in range(recs):
multdf.loc[geom,'geometry'] = row.geometry[geom]
outdf = outdf.append(multdf,ignore_index=True)
return outdf
#-------
#Explode the BritishColumbia MultiPolygon into its constituents
EBritishColumbia=explode_polygon(BritishColumbia)
#Loop over each individual polygon and get external coordinates
for index,row in EBritishColumbia.iterrows():
print 'working on polygon', index
mypolygon=[]
for pt in list(row['geometry'].exterior.coords):
print index,', ',pt
mypolygon.append(pt)
#See if any of the original grid points read from the netcdf file earlier
#lie within the exterior coordinates of this polygon
#pth.contains_points returns a boolean array (true/false), in the
#shape of 'points'
path=mpltPath.Path(mypolygon)
inside=path.contains_points(points)
#find the results in the array that were inside the polygon ('True')
#and set them to missing. First, must reshape the result of the search
#('points') so that it matches the mask & original data
#reshape the result to the main grid array
inside=np.array(inside).reshape(lon2d.shape)
i=np.where(inside == True)
mask[i]=1
print 'fininshed checking for points inside all polygons'
#mask now contains 0's for points that are not within British Columbia, and
#1's for points that are. FINALLY, use this to mask the original data
#(stored as 'fld')
i=np.where(mask == 0)
fld[i]=np.nan
#Done.