我处理包含气象数据(经度、纬度、降水、温度……)的虹膜立方体,并且我对计算特定区域(例如一个国家)的统计数据感兴趣。
这篇文章解释了如何用一个盒子(min lon、min lat、max lon、max lat)裁剪立方体,但我想更进一步,使用 shapefile 选择一个精确的区域。
这篇文章解释说可以使用与蒙版关联的 shapefile 来裁剪图像,但我不知道如何使它适用于我的虹膜立方体。
如果有人可以给我一个例子或解释我如何做到这一点,那将非常有用。
PS:我对python很陌生
使用例如Fiona读取 shapefile应该可以:
from shapely.geometry import MultiPoint
# Create a mask for the data
mask = np.ones(cube.shape, dtype=bool)
# Create a set of x,y points from the cube
x, y = np.meshgrid(cube.coord(axis='X').points, cube.coord(axis='Y').points)
lat_lon_points = np.vstack([x.flat, y.flat])
points = MultiPoint(lat_lon_points.T)
# Find all points within the region of interest (a Shapely geometry)
indices = [i for i, p in enumerate(points) if region.contains(p)]
mask[np.unravel_index(indices)] = False
# Then apply the mask
if isinstance(cube.data, np.ma.MaskedArray):
cube.data.mask &= mask
else:
cube.data = np.ma.masked_array(cube.data, mask)
这仅适用于 2D 立方体,但只需要针对更高维度进行调整,以便蒙版仅在纬度/经度维度上。
实际上,我最近在CIS中实现了这种行为,这样你就可以做cube.subset(shape=region)
这可能对你来说更容易。