Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

split_polygons_experimental failing on polygons crossing x or y datum #53

Open
thomas-fred opened this issue Jun 16, 2023 · 0 comments
Open
Labels
bug Something isn't working

Comments

@thomas-fred
Copy link

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

@thomas-fred thomas-fred added the bug Something isn't working label Jun 16, 2023
@thomas-fred thomas-fred changed the title split_polygon failing on polygons crossing x or y datum split_polygons_experimental failing on polygons crossing x or y datum Jun 19, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

No branches or pull requests

1 participant