Обрезать определенную область в кубе IRIS, используя шейп-файл

Question

Я работаю с кубиками радужной оболочки, содержащими метеорологические данные (долгота, широта, осадки, температура,...), и меня интересует расчет статистики в определенных областях (например, в стране).

В этом посте объясняется, как обрезать куб с помощью прямоугольника (min lon, min lat, max lon, max lat), но я бы хотел сделать шаг вперед и выбрать точную область, используя шейп-файл.

В этом посте объясняется, что можно обрезать изображение, используя шейп-файл, связанный с маской, но я не знаю, как заставить это работать для моих кубиков радужной оболочки.

Если бы кто-нибудь мог дать мне пример или объяснить, как это сделать, это было бы очень полезно.

PS: я довольно noobie с питоном

Answer 1

Прочитав шейп-файл, используя, например, Fiona, что-то вроде этого должно работать:

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-кубов, но просто требует настройки для более высоких размеров, чтобы маска была только над размерами широты и долготы.

Я фактически реализовал это поведение в СНГ недавно, чтобы вы могли сделать cube.subset(shape=region) что может быть проще для вас.