-
Notifications
You must be signed in to change notification settings - Fork 551
Description
The __eq__ function in CRS can be very slow if using CRSs are not from the EPSG authority, and in particular if they are custom CRSs as uncovered in cogeotiff/rio-tiler#747 (and before in davenquinn/mars-tiler#2)
The current implementation in
Lines 923 to 947 in 81ba490
| def __eq__(self, other): | |
| cdef OGRSpatialReferenceH osr_s = NULL | |
| cdef OGRSpatialReferenceH osr_o = NULL | |
| cdef CRS crs_o | |
| try: | |
| crs_o = CRS.from_user_input(other) | |
| except CRSError: | |
| return False | |
| epsg_s = self.to_epsg() | |
| epsg_o = crs_o.to_epsg() | |
| if epsg_s is not None and epsg_o is not None and epsg_s == epsg_o: | |
| return True | |
| else: | |
| try: | |
| osr_s = exc_wrap_pointer(OSRClone(self._osr)) | |
| osr_o = exc_wrap_pointer(OSRClone(crs_o._osr)) | |
| return bool(OSRIsSame(osr_s, osr_o) == 1) | |
| finally: | |
| _safe_osr_release(osr_s) | |
| _safe_osr_release(osr_o) |
From the discussion on rio-tiler I made some microbenchmarks pasted below that shows this slowdown as a symptom. What precisely needs to change (or should change) in CRS.__eq__ I am less sure of, but one potential work around would be to just compare the wkt property (which is cached) first if the other is a CRS class.
I am very happy to implement that change or other changes if needed.
Another possibility would be to update the conventions around to_epsg to allow the match to other authorities (so maybe a new function would be to_code or something).
I do have a PR in #3206 to add the IAU as a matching authority to to_authority which could be expanded to do that.
Steps to reproduce the problem.
In [57]: import rasterio as rio
In [58]: e = rio.CRS.from_epsg(4326)
In [59]: m = rio.CRS.from_user_input('IAU_2015:49900')
In [60]: m
Out[60]: CRS.from_wkt('GEOGCS["Mars (2015) - Sphere / Ocentric",DATUM["Mars (2015) - Sphere",SPHEROID["Mars (2015) - Sphere",3396190,0,AUTHORITY["IAU","49900"]],AUTHORITY["IAU","49900"]],PRIMEM["Reference Meridian",0,AUTHORITY["IAU","49900"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["IAU","49900"]]')
In [61]: %timeit e == e
77.6 ns ± 11.4 ns per loop (mean ± std. dev. of 7 runs, 10,000,000 loops each)
In [62]: %timeit m == m
39.7 μs ± 6.81 μs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
In [63]: %timeit m.to_wkt() == m.to_wkt()
7.95 μs ± 944 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
In [64]: %timeit m.wkt == m.wkt
36.4 ns ± 6.9 ns per loop (mean ± std. dev. of 7 runs, 10,000,000 loops each)
In [72]: m2 = rio.CRS.from_wkt('GEOGCS["Mars (2015) - Sphere XY",DATUM["Mars (2015) - Sphere",SPHEROID["Mars (2015) - Sphere",3396190,0]],PRIMEM["Reference Meridian",0]
...: ,UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Longitude",EAST],AXIS["Latitude",NORTH]]')
In [73]: m2
Out[73]: CRS.from_wkt('GEOGCS["Mars (2015) - Sphere XY",DATUM["Mars (2015) - Sphere",SPHEROID["Mars (2015) - Sphere",3396190,0]],PRIMEM["Reference Meridian",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Longitude",EAST],AXIS["Latitude",NORTH]]')
In [74]: %timeit m2 == m2
20.8 ms ± 23.3 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
In [75]: %timeit m2.wkt == m2.wkt
36.3 ns ± 6.81 ns per loop (mean ± std. dev. of 7 runs, 10,000,000 loops each)Environment Information
rasterio info:
rasterio: 1.4.1
GDAL: 3.9.2
PROJ: 9.5.0
GEOS: 3.13.0
PROJ DATA: /home/aannex/micromamba/envs/marspre/share/proj
GDAL DATA: /home/aannex/micromamba/envs/marspre/share/gdal
System:
python: 3.12.7 | packaged by conda-forge | (main, Oct 4 2024, 16:05:46) [GCC 13.3.0]
executable: /home/aannex/micromamba/envs/marspre/bin/python3.12
machine: Linux-6.5.0-1027-oem-x86_64-with-glibc2.35
Python deps:
affine: 2.4.0
attrs: 24.2.0
certifi: 2024.08.30
click: 8.1.7
cligj: 0.7.2
cython: None
numpy: 2.1.2
click-plugins: None
setuptools: 75.1.0
Installation Method
For example: from conda-forge via micromamba