région définie découper un cube dans IRIS en utilisant un shapefile

Question

Je manipule des cubes d'iris contenant des données météorologiques (longitude, latitude, précipitations, température, ...) et je suis intéressé par le calcul des statistiques dans les zones définies (par exemple un pays).

Cette post explique comment recadrer le cube avec une boîte (min lon, min lat, max lon, max lat) mais je voudrais aller plus loin et sélectionner une zone précise en utilisant un fichier de formes.

Cette post explique qu'il est possible de rogner une image en utilisant un shapefile associé à un masque, mais je ne sais pas comment je peux le faire fonctionner pour mes cubes d'iris.

Si quelqu'un pouvait me donner un exemple ou me expliquer comment faire ce serait très utile.

PS: Je suis tout à fait noobie avec python

Answer 1

Après avoir lu le shapefile en utilisant par exemple Fiona quelque chose comme cela devrait fonctionner:

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) 

Cela ne fonctionne que pour les cubes en 2D, mais juste pour les dimensions doit être modifié de façon plus que le masque est seulement sur les latitude/longitude dimensions.

J'ai récemment implémenté ce comportement dans CIS récemment afin que vous puissiez faire cube.subset(shape=region) qui pourrait être plus facile pour vous.