nismod / snail

spatial networks impact assessment library 🐌
https://nismod.github.io/snail/
MIT License
9 stars 1 forks source link

split_polygons_experimental failing on polygons crossing x or y datum #53

Closed thomas-fred closed 4 months ago

thomas-fred commented 1 year ago

Trying to split polygons which straddle x=0 or y=0 by a raster grid with snail.core.intersection.split_polygon_experimental, sometimes triggers RuntimeError from here.

Code to reproduce for 0.3.1 is something like:

#!/usr/bin/env python
# coding: utf-8

"""
Failing test case for snail==0.3.1

Try to split a polygon by a raster grid, where that polygon crosses a coordinate zero line.

Fails with `RuntimeError: Expected even number of crossings on gridline.`
"""

import geopandas as gpd
import shapely

import snail
from snail.intersection import (
    Transform,
    prepare_polygons,
    split_polygons_experimental
)

bad_poly = shapely.wkt.loads(
    (
        'POLYGON (('
        '-0.0062485600499826 51.61041647955, '
        '-0.0062485600499826 51.602083146149994, '
        '0.0020847733500204 51.602083146149994, '
        '0.0020847733500204 51.61041647955, '
        '-0.0062485600499826 51.61041647955'
        '))'
    )
)
bad_gdf = gpd.GeoDataFrame(geometry=[bad_poly])

# mimicking cli.split
transform = Transform(
    crs=(
        'PROJCS["World_Mollweide",GEOGCS["WGS 84",DATUM["WGS_1984",'
        'SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],'
        'AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],'
        'UNIT["Degree",0.0174532925199433]],PROJECTION["Mollweide"],'
        'PARAMETER["central_meridian",0],PARAMETER["false_easting",0],'
        'PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],'
        'AXIS["Easting",EAST],AXIS["Northing",NORTH]]'
    ),
    width=36082,
    height=18000,
    transform=(1000.0, 0.0, -18041000.0, 0.0, -1000.0, 9000000.0)
)
prepared = prepare_polygons(bad_gdf)

# will raise RuntimeError; number of gridline crossings is not even
splits = split_polygons_experimental(prepared, transform)

See below for an example with gridfinder 'targets' and a population raster.

snail_meridian_failures_context

Not all polygons overlapping Greenwich meridian fail to split, but many do. snail_meridian_failures_closeup