Any way to init a projection directly from PROJCS strings?

[dannguyen]

I know this isn't StackOverflow but just wanted to make sure I wasn't missing something obvious that I didn't see while reading through the source.

I have a PROJCS string that looks like this:

 'PROJCS["NAD_1983_StatePlane_New_York_Long_Island_FIPS_3104_Feet",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",984250.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-74.0],PARAMETER["Standard_Parallel_1",40.66666666666666],PARAMETER["Standard_Parallel_2",41.03333333333333],PARAMETER["Latitude_Of_Origin",40.16666666666666],UNIT["Foot_US",0.3048006096012192]]'

I've seen other projection libraries (such as proj4 for node) handle this, and so I was wondering if pyproj.Proj() could as well? Trying to construct it from that raw string gives me this error:

          RuntimeError: b'projection not named'

Sorry for the support question. I don't know much about GIS systems...pyproj.Proj() has worked just fine for me doing this:

          pyproj.Proj(init="EPSG:4269")

and:

         pyproj.Proj("""
            +proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333 +lat_0=40.16666666666666 +lon_0=-74
            +x_0=984249.9999999999 
            +y_0=0 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs""")

And so reading from PROJCS type strings (or is it referred to as WKT format?) seems like one more translation. But I don't know much about GIS beyond knowing that I'm trying to convert between two systems, so I apologize if I'm asking about something ridiculously apples-to-oranges.

[dannguyen]

OK, I actually solved my original problem. I didn't realize pyproj.Proj() explicitly ignores the +to_meter and +units=us-ft parameters:

import pyproj
x, y = (984533, 213998)
projstr = """+proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333 +lat_0=40.16666666666666 +lon_0=-74 
+x_0=300000 +y_0=0 +ellps=GRS80 +datum=NAD83 
+to_meter=0.3048006096012192 +units=us-ft +no_defs"""

myproj = pyproj.Proj(projstr, preserve_units=True)
myproj(x, y, inverse=True)
# (-73.99897854592253, 40.75405155302713)

badproj = pyproj.Proj(projstr)
badproj(x, y, inverse=True)
# (-65.75200831429146, 41.80318458470038)

I imagine that it's not possible to make preserve_units=True by default as to not break previous behavior? But is it worth putting a note in the documentation so that this is clearer...particularly for casual GIS users (like me :) ) who assume that copying-pasting a the proj4 string from places like spatialreference is all that's needed?

This was alluded to in this issue: #40

[micahcochran]

You aren't the only one that was tripped up by that flag. I don't really like the default behavior of the preserve_units flag, but this is in the API. If it were changes, it would probably cause bugs in other code that doesn't set the flag to False.

The current code defaults to adding a +units=m parameter even when it doesn't make sense, for instance:

In [7]: from pyproj import Proj

In [8]: p = Proj(init='epsg:4326')

In [9]: p.srs
Out[9]: '+units=m +init=epsg:4326 '

In [10]: p.is_latlong()
Out[10]: True

EPSG:4326 is aka WGS84, which is usually units decimal degrees--I think pyproj has some radian conversions that it uses. It seems that PROJ.4 ignores +units=m when it doesn't make sense.

I assume a pull request would be welcome that better documents the preserve_units. The documentation is generated directly from the code, so edit the init.py file.

[dannguyen]

I think the function documentation works for what it is, I mean in terms of following a standard docstring. What the problem for me -- and I imagine most new users -- is that Googling around for examples of pyproj yields straightforward examples that don't mention preserve_units...In retrospect, most of these examples were from authors outside the United States, i.e. people who don't have to deal with America's/England's ft/meters crapfest :)

I was thinking of having something on either the README.md or a separate documentation page full of examples. I'd be happy to write a few...maybe throw it in the Wiki, which can then be later converted to a proper documentation page? Just pointing out that while preserve_units is documented in the code...and it is fairly straightforward to grok once you know of its existence...virtually none of the Googlable examples mention it, and for some people (again I'm new to GIS so I may just have an idiotic way of thinking about things), the instinct might be to check for errors in the proj4 string itself, rather than suspecting pyproj.Proj(). Having an upfront example that mentions it would cut that confusion time significantly.

[micahcochran]

Still having three different measurements for feet is great, which is clearly a distinct advantage for feet over meters. 🙈 Feet only has one name and meters can also be spelled, metres, which my browser correctly flags as the incorrect spelling 🙉 . End of sarcasm.

How about this. To view, either download the HTML or download the ipynb. For the ipynb, you'll need ipython (with notebooks) and folium installed to see it.

Hopefully with the HTML, I don't need to upload anything other files.

Also, the x/y order versus lat/long tends to confuse some people.

[dave-nm]

@micahcochran Could we do a 2.0 release that breaks the API and changes preserve_units to default to True?

[micahcochran]

I'm hesitant about doing that. The Proj.init appends a +units=m when preserve_units is False, even if you are using +proj=latlong (proj4 ignores that). There might be a good reason for doing this, if you are actually crafting proj4 strings by hand--so to speak. People using feet are usually not scientists.

pyproj goes back to atleast 2006, which it added support for Python 2.5 and supported earlier Python version! Sure, the library has a few design decisions that might not be considered Pythonic by today's standards. It'd be nice if the transform function would take a list of coordinate tuples. I'd prefer not to break API capability, if possible.

My hesitance is that it could cause some old code to break. How would you notify the user that preserve_units defaults to True? You are not deprecating behavior that is not longer considered the correct way to do things.

In summary, I'm hesitant to make that kind of change.

[dannguyen]

I'm a little more experienced with GIS now, at least from when I first filed this report :)

To give you some context, I was trying to do geospatial analysis of the NYPD Stop and Frisk data found here:

http://www.nyc.gov/html/nypd/html/analysis_and_planning/stop_question_and_frisk_report.shtml

Here's an example of what such analysis used for, journalistically: http://johnkeefe.net/nypd-stop-frisk-data-for-you

The NYPD, and all other NYC agencies, use the New York-Long Island State Plane Coordinate System. I'm not a ArcGIS or QGIS user, so I don't know what the typical convention is to find that info, but as a casual GISer, I just google it and eventually find my way to spatialreference.org:

http://www.spatialreference.org/ref/esri/102718/

The "Proj4" link produces this text string, which is what I referenced in my original issue report:

      +proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs 

So it wasn't a string that I was hand-encoding, but what I had thought the normal convention for GIS to be: give a string to the projection formula to declare how I want coordinates to be translated.

If I instead invoke Proj using the name of the projection, e.g. ESRI:102718, I get the same unexpected mistranslation:

F = pyproj.Proj(init="ESRI:102718")
T = pyproj.Proj(init="ESRI:102718", preserve_units=True)
x = 1056308
y = 182711

>>> F(x, y, inverse=True)
(-64.93197150112687, 41.459120016286946)
>>> T(x, y, inverse=True)
(-73.74025103331826, 40.66788294565469)

I don't know enough about the rest of the GIS scene to know how other projections work, but I can say that for many American encoding systems, feet seems to be the standard, and it's not really about what's best scientific practice...I think we can all agree that the metric system is the way to go :)

I think it's fine for pyproj to have the opinion that projections should be using meters as units, as a default. But my quibble is that even when the user does the explicit, painful thing, which is to declare the units in the string, Proj will override that explicit declaration.

That said, it's an open question whether it is worth a breaking change to fix the default preserve_units parameter. If there is going to be a 2.0 release that makes other breaking changes, then sure, maybe it's worth it. But is it worth making a major point release just for this? Again, I'm just a casual user so I have no idea whether that makes sense, given other people's use case.

Maybe one mitigation strategy could be to have Proj() check if units is mentioned in the init definition. If so, then issue a warning that Proj() will override it with +units=m, and that if the user doesn't like that, then they should use preserve_units=True. With that, and a clarification in the documentation, that would seem to effectively mitigate errors for new users without breaking existing scripts.

[vsasseville]

Hello, I would be interested in the answer of your original question. How can you extract the projection from a WKT string such as:
WKT_string = "PROJCS["NAD_1983_Terr4M",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-105.0],PARAMETER["Standard_Parallel_1",60.0],PARAMETER["Standard_Parallel_2",75.0],PARAMETER["Latitude_Of_Origin",60.0],UNIT["Meter",1.0]]"

And have it in a pyproj object like pyproj.Proj(wkt_string)

I can't find how to do this anywhere and it looks like it should be simple. Thx

For broader context: I want to export the projection from a shapefile (or .prj file) in a python object to use pyproj.transform with it.

[micahcochran]

@vsasseville
PyCRS does a decent job converting Well Know Text to a PROJ.4 string. It is all python.

GDAL/OGR Python's SpatialReference functions can convert between WKT and PROJ.4 strings (and a few more), but you have to have GDAL installed and that instance of Python has to be able to get to it.

[vsasseville]

Didn't know about PyCRS, thanks a lot!

[snowman2]

WKT projection strings will be supported in the pyproj 2.0.0 release. See: #152

[jswhit]

Master has been updated to 2.0.0, which requires proj4 6.0.0, and changes the default value of preserve_units (thanks to @snowman2).

[snowman2]

2.1.0 is released and should resolve this issue. Thanks for the report!