0

I handle iris cubes containing meteorological data (lon, lat, precipitation, temperature,...) and I am interested in calculating statistics in defined areas (for example a country).

This post explains how to crop the cube with a box (min lon, min lat, max lon, max lat) but I would like to go a step further and select a precise area using a shapefile.

This post explains that it is possible to crop an image using a shapefile associated to a mask, but I don't know how I can make it work for my iris cubes.

If somebody could give me an example or explain me how to do that it would be very useful.

PS: I am quite noobie with python

Community
  • 1
  • 1
Zoran
  • 309
  • 4
  • 9

1 Answers1

1

Having read the shapefile using e.g. Fiona something like this should work:

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)

This only works for 2D cubes, but just needs tweaking for higher dimensions so that the mask is only over the lat/lon dimensions.

I actually implemented this behaviour in CIS recently so that you can do cube.subset(shape=region) which might be easier for you.

Duncan WP
  • 830
  • 5
  • 19