Comments (13)
@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.
from pyproj.
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
from pyproj.
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.
from pyproj.
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.
from pyproj.
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.
from pyproj.
@micahcochran Could we do a 2.0 release that breaks the API and changes preserve_units
to default to True
?
from pyproj.
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.
from pyproj.
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.
from pyproj.
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.
from pyproj.
Didn't know about PyCRS, thanks a lot!
from pyproj.
WKT projection strings will be supported in the pyproj 2.0.0 release. See: #152
from pyproj.
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).
from pyproj.
2.1.0 is released and should resolve this issue. Thanks for the report!
from pyproj.
Related Issues (20)
- DeprecationWarning: pkg_resources is deprecated as an API. HOT 1
- Support for deprecated and non-deprecated objects HOT 3
- Invert lat lon "Transformer" HOT 1
- Internal Proj Error: proj_create: unrecognized format / unknown name HOT 4
- Inverted x/y for transformation from WSG84 to LV95 (EPSG:2056) HOT 3
- muslinux wheel missing HOT 1
- test_azimuthal_equidistant fails with PROJ 9.4.0 HOT 1
- test failure in test_crs_cf.py HOT 2
- Add proj_query_geodetic_crs_from_datum HOT 2
- 3.6.1: pytest fails in `test/crs/test_crs_cf.py::test_azimuthal_equidistant` unit HOT 1
- A question about pyproj/_datadir.pyx HOT 2
- geometry_area_perimeter limitations not documented HOT 1
- Calculate area of arbitrary polygon on Earth's surface HOT 4
- Interrupted Mollweide gives wrong result on inverse projections of some lat/lons. HOT 1
- test_datum__from_epsg fails with PROJ 9.4.1 HOT 4
- Proj4 export from a bound CRS object fails HOT 4
- Antipodes don't match when calculated numerically HOT 7
- Bad_alloc with pyroj import in QGIS plugin
- Conversion of array to scalar when transforming 1 element arrays leads to DeprecationWarning HOT 3
- Small memory leak HOT 3
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from pyproj.