Pyproj test for two equal projections

[fmaussion]

pyproj doesn't seem to check if projections are equal when asked to transform:

x = np.random.randn(1e7) * 60
y = np.random.randn(1e7) * 60
start_time = time.time()
xx, yy = pyproj.transform(pyproj.Proj(init='epsg:26915'), pyproj.Proj(init='epsg:26915'), x, y)
print(end_time-time.time())

Takes about 3 seconds on my laptop. I wonder if it would be possible for pyproj to catch these?

Would the following workaround work or is it dangerous?

def my_transform(p1, p2, x, y):
    if p1.srs == p2.srs:
        return x.copy(), y.copy() # or no copy if kwd nocopy=True
    return pyproj.transform(p1, p2, x, y)

Thanks,

Fabien

[micahcochran]

I like the idea of a few lines of code detecting if the projections are the same. I am unsure if the built-in transform() should do that checking.

Your work around will work for the above case because you defined both projections in exactly the same way. It is not a bad way to test. Look at my example using the same projection defined in two different ways.

p1 = pyproj.Proj(init='epsg:26915')
p2 = pyproj.Proj('+proj=utm +zone=15 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0')
if p1.srs == p2.srs:
  print("They are equal???")     # this shouldn't happen
else:
  print("Test fails.")

In order to make a projection equivalence test, PyProj would either have to expose pj_get_def() from Proj.4 or create a new function Proj.__eq__() that to check if the projections are equivalent and most likely use the `pr_get_def() .

I like the idea of a Proj.__eq__() so you can test p1 == p2. Then, your code could test == and skip the transform. There other useful application for that projection equivalency of function.

[fmaussion]

Yes I see. I would find a Proj.__eq__() very useful in general, yes! How complicated would it be to implement?

[micahcochran]

Pseudocode of how I would make a Proj.__eq__().

1. Run pj_get_def() on a projection
2. convert Proj4 string to Fiona.crs dictionary representation, from_string()
3. remove those key from dictionary that are for metadata, such as the 'init' key
4. Do steps 1-3 on a 2nd projection.
5. compare parameters that a functionally equivalent or convert them, such as '+to_meters=.3048' and '+units=ft'
6. return Boolean value of p1_dict == p2_dict

This seems fairly easy to me.

Below are a few more notes:
Running pj_get_def() for the projection string '+init=epsg:26915'
returns the string ' +init=epsg:26915 +proj=utm +zone=15 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0'

Fiona has crs.py, which has the from_string(), which converts a above Proj.4 strings into dictionaries. If the licenses will allow, I would just copy from_string() into pyproj.

I don't know Cython, which pyproj uses to interface with the Proj.4, but Cython seems fairly straightforward.

[micahcochran]

Just to add one more thing. The goal of this code is very similar to the already existing OGR function IsSame(), (C++ code) which can be accessed in Python through osgeo.osr library.

If you need that type of functionality right now, I'd use that library. If you are already using Fiona, using osr shouldn't be any more of a burden to your users. PyProj probably should not be relying upon OGR.

[fmaussion]

@micahcochran sorry it took me so long to get back to you. Indeed your suggestion is working perfectly (see test case below). I agree with you that pyproj should avoid depending upon OGR (the fact that pyproj is lightweight is really great). I just changed to the title and keep the issue open for later use. A "homemade" p1.is_same(p2) would still be very useful.

EDIT: or even more pythonic: p1.__eq__(other).

IsSame using OGR:

import pyproj
from osgeo import osr 

# Pyproj version
p1 = pyproj.Proj(init='epsg:26915')
p2 = pyproj.Proj('+proj=utm +zone=15 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0')
if p1.srs == p2.srs:
      print("They are equal")
else:
      print("Test fails.")  # This is what happens

# OGR version
s1 = osr.SpatialReference()
s1.ImportFromProj4(p1.srs)
s2 = osr.SpatialReference()
s2.ImportFromProj4(p2.srs)
if s1.IsSame(s2):
      print("They are equal")   # This is what happens
else:
      print("Test fails.")

[micahcochran]

No problem @fmaussion

In my limited additional testing and OGR/GDAL seems to fail sometimes, but that does work better for more cases than p1.srs == p2.srs.

This would be extending pyproj beyond the capabilities of proj.4's public facing API, which pyproj mirrors fairly closely. OGR/GDAL depends upon proj.4 for its coordinate transformations, but extends its capabilities with the IsSame and other functions.

I'm going to make a suggestion that's not quite going to solve your problem. PyCRS is a python library has some of the workings to understand projection parameters. It is already a useful library to convert between WKT and proj4 projection representations. It seems like a projection equality comparison could compliment it. It might not be too hard extend it. Note: PyCRS still has a few rough edges.

[fmaussion]

Thanks a lot @micahcochran ! I'll have a look at PyCRS, although I think that the simple solution above will largely fulfill my needs.

I'll close this issue then.