Skip to content

reproject: spherical-ellipsoid CRSes (MODIS sinusoidal) take the WGS84 fast path, ~19-24 km error #3275

Description

@brendancol

Describe the bug

The numba/CUDA fast-path dispatch in xrspatial/reproject/_projections.py applies WGS84 ellipsoid constants to CRSes that are not on a WGS84-compatible ellipsoid. Two holes:

  1. _is_wgs84_compatible_ellipsoid only looks at the ellps and datum keys of crs.to_dict(), treating empty values as WGS84. A spherical CRS defined by explicit axes is normalized by PROJ to an R key (no ellps), so it passes the check. The standard MODIS sinusoidal grid is exactly this case:

    import pyproj
    from xrspatial.reproject._projections import try_numba_transform
    modis = pyproj.CRS("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs")
    # try_numba_transform(crs4326, modis, ...) takes the fast path

    The WGS84 meridional-arc kernel then runs against a sphere. Measured against pyproj: 18.9 km max error at mid latitudes. Any reproject() between EPSG:4326 and a MODIS sinusoidal grid silently returns data shifted by that much. The CUDA dispatch (try_cuda_transform) shares the helper, so the GPU path is wrong the same way.

  2. _aea_params and _cea_params have no ellipsoid check at all (every other param extractor calls _is_wgs84_compatible_ellipsoid). A spherical AEA or CEA definition runs through the WGS84 kernels:

    +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +a=6370997 +b=6370997   -> up to 23.8 km error
    +proj=cea +lat_ts=30 +a=6371228 +b=6371228                                      -> up to 23.8 km error
    

    These slip past the datum-shift guard because ellps=sphere is not in _DATUM_PARAMS.

Expected behavior

CRS pairs the kernels cannot represent should return None from try_numba_transform / transform_points / try_cuda_transform and fall back to pyproj, the same way OSGB36/NAD27 pairs already do (#2651).

Additional context

Found by /sweep-accuracy against the reproject module. The fix needs _is_wgs84_compatible_ellipsoid to reject R-defined spheres and explicit a/b axes that differ from WGS84, and _aea_params/_cea_params to call the guard. WGS84/GRS80 pairs must keep the fast path (GRS80 axes differ from WGS84 by ~1e-4 m, so an exact-match check is too strict).

Activity

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    bugSomething isn't workinggpuCuPy / CUDA GPU supportsweep-accuracyFound by /sweep-accuracy

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions