Bugfix shape masking

[hsteptoe]

🚀 Pull Request

Description

Fixes #6126 by exposing the shape geometry to the user. Adds to docs and docstrings to make it clearer how to use this.

---

Consult Iris pull request check list

---

Add any of the below labels to trigger actions on this PR:

• benchmark_this Request that this pull request be benchmarked to check if it introduces performance shifts

[pp-mo]

From @SciTools/peloton : discussed what we are seeing here and in #6126.
We think this is looking useful, but the two of you seem to have it covered @hsteptoe @acchamber ?
Please shout if you need help.

[hsteptoe]

Discussion with @acchamber notes:

• Don't move imports out of functions as this creates circular import errors
• Check target_crs assignment for occasion when cube doesn't have a CRS
• Re-write conversion of 0-360 -> -180-180 coords in case of degree based CRSs

[hsteptoe]

What are @SciTools/peloton (@pp-mo) thoughts about adding rasterio as a Iris dependency?

This would facilitate better handling of shapes to rasters (which is essentially the problem we're solving) with the added bonus that we could also mask to other shape types, like lines, to extract a trajectory from a Cube, for example.

[acchamber]

I've had a few other discussions with other users now where they can't use this function with OSGB shapefiles which this PR would fix. Where are we on being able to include it in the next release?

@hsteptoe @pp-mo

[hsteptoe]

I've had a few other discussions with other users now where they can't use this function with OSGB shapefiles which this PR would fix. Where are we on being able to include it in the next release?

@hsteptoe @pp-mo

I was planning on going to the AVD surgery this week to discuss some changes and/or make the case for rasterio

[trexfeathers]

Thanks for your patience everyone. It's sometimes a struggle to balance everything and we just haven't had an opportunity to discuss this. Coming to the Surgery (UK Met Office) is an ideal next step.

[stephenworsley]

@pp-mo, @trexfeathers and @stephenworsley have discussed this at the Surgery and agree in principle that we could consider adding rasterio as an optional dependency.

[trexfeathers]

From @SciTools/peloton: would @hsteptoe and @acchamber appreciate any more input from core devs at this point?

[hsteptoe]

From @SciTools/peloton: would @hsteptoe and @acchamber appreciate any more input from core devs at this point?

I just need to find more time, which is in short supply as we get near the end of FY...

[hsteptoe]

This is now ready for some initial (alpha?) testing 🎉. @acchamber are you available to try out the new features? (Updates to docs etc. will follow if this is successful) Tests are incomplete are should not be review yet.

This is a fairly substantial re-write. New/improved features include:

• Automatic CRS transforms. Will handle situation where the cube and geometry have different CRSs by transforming the geometry.
• Handles all shape types, including Points, Lines and MultiPolygons.

There are some breaking changes, principally moving away from having a minimum_weight keyword, replacing it with an all_touched keyword. The rational here is for consistency across shape types. Minimum weight doesn't make sense for Points and Lines. Instead, Line pixilation is based on Bresenham's line algorithm. Polygon pixilation is based on either: (a) pixel centres that intersect the polygon, or (b) all pixels touched by the polygon irrespective of the area of intersection.

My informal tests so far show working behaviour for iris test data from rotated_pole.nc, toa_brightness_stereographic.nc, E1_north_america.nc, air_temp.pp. Know limitation (not-working) for semi-structured model grids, such as orca2_votemper.nc. Some test shapefiles 👉test_shapefiles.zip

[hsteptoe]

@trexfeathers @bjlittle I think this is ready for a first review... in particular I would appreciate reviews/help with

• Checking the logic of the mutually exclusive kwargs minimum_weight and all_touched. Not sure if I've done this optimally
• Explaining some of the CI/pre-commit failures... the numpydoc errors in particular are a bit opaque to me...
• Checking the best place to set a crucial PROJ env variable (see b52a0d7)

[mo-DanCubbon]

@trexfeathers @bjlittle I think this is ready for a first review... in particular I would appreciate reviews/help with

• Checking the logic of the mutually exclusive kwargs minimum_weight and all_touched. Not sure if I've done this optimally
• Explaining some of the CI/pre-commit failures... the numpydoc errors in particular are a bit opaque to me...
• Checking the best place to set a crucial PROJ env variable (see b52a0d7)

re. the PROJ variable, the last time I saw this come up was <URL> . Instead of having an SSL certificate problem, the proj transformation "succeeds" but actually fails silently and returns infinities.

[trexfeathers]

@trexfeathers @bjlittle I think this is ready for a first review... in particular I would appreciate reviews/help with

• Checking the logic of the mutually exclusive kwargs minimum_weight and all_touched. Not sure if I've done this optimally
• Explaining some of the CI/pre-commit failures... the numpydoc errors in particular are a bit opaque to me...
• Checking the best place to set a crucial PROJ env variable (see b52a0d7)

re. the PROJ variable, the last time I saw this come up was <URL> . Instead of having an SSL certificate problem, the proj transformation "succeeds" but actually fails silently and returns infinities.

@mo-DanCubbon only two-thirds of Iris devs will be able to view this link since it's internal to the UK Met Office. If there is any more detail you'd like to be visible it'll need to be posted verbatim.

[mo-DanCubbon]

@trexfeathers @bjlittle I think this is ready for a first review... in particular I would appreciate reviews/help with

• Checking the logic of the mutually exclusive kwargs minimum_weight and all_touched. Not sure if I've done this optimally
• Explaining some of the CI/pre-commit failures... the numpydoc errors in particular are a bit opaque to me...
• Checking the best place to set a crucial PROJ env variable (see b52a0d7)

re. the PROJ variable, the last time I saw this come up was <URL> . Instead of having an SSL certificate problem, the proj transformation "succeeds" but actually fails silently and returns infinities.

@mo-DanCubbon only two-thirds of Iris devs will be able to view this link since it's internal to the UK Met Office. If there is any more detail you'd like to be visible it'll need to be posted verbatim.

Details in that link were not super important, only that the network had to be switched off. We've investigated this a bit more, and have turned our PROJ_NETWORK back to on... but now we're getting the right answer. I'm assuming Proj has cached the necessary grid files to perform the transformation correctly now. I haven't since been able to reproduce the issue unfortunately but you could try with this:

from pyproj import Transformer
from cartopy import crs as ccrs
import pyproj.network

EASTING_NORTHING_EXAMPLE = (300000, 100000)

if __name__ == "__main__":
    print(f"Is network enabled: {pyproj.network.is_network_enabled()}")
    # First transform to OSGB
    osgb_crs = ccrs.OSGB()
    wgs84_crs = ccrs.PlateCarree()
    transformer = Transformer.from_crs(osgb_crs, wgs84_crs, always_xy=True)
    latlon = transformer.transform(*EASTING_NORTHING_EXAMPLE)

    print(f"coordinates: {latlon}")

    # Now disable network
    pyproj.network.set_network_enabled(active=False)
    print(f"Is network enabled after disabling: {pyproj.network.is_network_enabled()}")

    # instantiate the transformer again
    transformer = Transformer.from_crs(osgb_crs, wgs84_crs, always_xy=True)
    latlon_no_network = transformer.transform(*EASTING_NORTHING_EXAMPLE)

    print(f"coordinates without network: {latlon_no_network}")

[mo-DanCubbon]

Just a question about the terminology used in this area. One of the functions being modified is mask_cube_from_shapefile is really masking from a geometry (which can contain multiple disconnected shapes) and not a shapefile (which could contain multiple geometries). I would guess it would be difficult to change the name of the function, but perhaps it is worth clarifying this in the documentation as it was something that initially confused me.

[hsteptoe]

Yes, that's a fair comment. The primary consideration here was not to introduce breaking changes by changing the function name.

I would be in favour of renaming for consistency with the shapely library, which underpins the geometry handling, and either go for mask_cube_from_shape() or mask_cube_from_geometry() which would recognise that iris is agnostic to how geometries are loaded in the first place (they could also come from geoJSON or .kml files for example), but @trexfeathers would need to approve the change... could we utilise iris.FUTURE here?

[trexfeathers]

This has been raised before, but during design @acchamber assured me that it would be much more meaningful to our user base if we used the word shapefile. Would you entertain the possibility that those of us discussing this on GitHub are more technical than the average user?

[hsteptoe]

Yes, this degree of technical pedantry may be unnecessary...but arguably it provides clarity. The simplest solution is just to add some clarification to the docs (although IMO the typing already makes it clear what the function expects). I suspect that the current convention is only meaningful to those folk coming from ASCEND (which limits it to internal MO), as ASCEND did also handle shapefile loading, which iris does not.

[acchamber]

I can't remember previous discussuions on this, so I may be contractdicting past self, but I don't see a big reason not to to use mask_cube_from_shape() instead - @hsteptoe is right in that it is not the shapefile itself being used, but a shapely geometry object. I think mask_cube_from_shape() makes it purpose clear to scanners of the docs while also being more correct - but that probably does require a depreciation warning for a switch in a few versions time, similar to how numpy handled this.

[trexfeathers]

OK well if everyone is cool with this, it's not as bad as you might think. You can introduce the new signature right now, and the outdated signature just calls the new one, with a deprecation warning as you say.

We are talking more about Iris v4.0 these days so the warnings won't be around forever!

[hsteptoe]

I think my preference would be mask_cube_from_shape()... @mo-DanCubbon @acchamber thoughts?

[mo-DanCubbon]

I think that would be the most understandable yes

[acchamber]

agreed