# Tag Info

5

If you're using pyproj2, it's much easier to use a Transformer. Here's an example: import pyproj from shapely.ops import transform project = pyproj.Transformer.from_proj( pyproj.Proj(init='epsg:4326'), # source coordinate system pyproj.Proj(init='epsg:26913')) # destination coordinate system # g1 is a shapley Polygon g2 = transform(project....

4

This is a great question and I will do my best to answer. To begin, the init style syntax is deprecated (https://pyproj4.github.io/pyproj/stable/gotchas.html#init-auth-auth-code-should-be-replaced-with-auth-auth-code). So, instead of CRS(init="epsg:4544"), you should use CRS("epsg:4544"). I discovered that sometimes when intialising a pyproj Proj ...

4

First, let's see what the mapping between a Y-X plane and a Lon-Lat plane would look like: Notes: Y-X-H coordinates are Cartesian ones. Since you named it Gaussian coordinates, it is suposed to Y and X are Transverse Mercator projections of the geodetic coordinates. I am assuming that Y coordinates are Eastings and X coordinates are Northings. ...

3

I expected this code to reproject the WGS84 point to a point on the NDFD grid, but it doesn't. Why not? Your code does reproject the WGS84 point to a point on the NDFD grid. Your problem is that rasterio is not picking up your GRIB file georeferencing and just returning the bounds as pixel coordinates from 0-ncols, 0-nrows (in fact your code should not ...

3

You could use an Azimuthal Equidistant projection centred on your origin point. In PROJ4 syntax, "+proj=aeqd +lat_0=1.234 +lon_0=5.678". That long/lat coordinate should project to (0,0) in that coordinate system, and points away from there should be in metres. You might want to add "+datum=WGS84 +units=m" but they may be default. Not sure what ellipsoid is ...

3

You can directly use shapely or GeoPandas but with 9888562 records It will take a long time to do (if you want a Progress bar during the pandas operations, you can use tqdm: ) 1) With your solution and the first 4 points import pandas as pd df = pd.DataFrame({'LAT':[47.9767,47.9803,47.9801,47.9798], 'LON':[-122.2450,-122.2458,-122.2472,-122.2465]}) ...

3

Use pyproj.CRS from pyproj import CRS CGCS2000metre = CRS('epsg:4544') print(CGCS2000metre.to_epsg()) # prints "4544" To get a CRS in the way you're already using pyproj, you can get a CRS object like this: from pyproj import Proj CGCS2000metre = Proj(init='EPSG:4544').crs print(CGCS2000metre.to_epsg()) # prints None print(CGCS2000metre.to_epsg(...

3

You need to set the environment variables to run the QGIS python outside. In this case my QGIS installation is in D:\QGIS, so you should only change OSGEO4W_ROOT. @ECHO OFF set OSGEO4W_ROOT=D:\QGIS call "%OSGEO4W_ROOT%\bin\o4w_env.bat" call "%OSGEO4W_ROOT%\bin\qt5_env.bat" call "%OSGEO4W_ROOT%\bin\py3_env.bat" path %OSGEO4W_ROOT%\apps\qgis\bin;%PATH% ...

2

Andrea Aime answered it in the geoserver-users list, which I quote here for the records: the syntax looks like a WKT2 representation, which GeoServer is not able to parse (code contributions or funding to get it going more than welcomed). I believe you'll have to rewrite it using WKT compatible sytanx instead (I'm not aware of an automatic ...

2

The is no solution to your problem, unfortunately. On one hand, as soon as you project coordinates on a 2 dimensional map, you need to choose between an equal-area projection (preserving the area) or a conformal projection (preserving local angles), but no projection has both. Of course, some projections make a compromise so that neither are totally ...

2

My guess is that your pyproj is too new. I reinstalled my setup on a different computer and had to roll back all the way to pyproj 1.9.6. I haven't dug deeper into it, but i think that pyproj has renamed something, as it is the line: from pyproj import Proj, transform, set_datapath that throws an exception. After changing to pyproj 1.9.6 it works. In ...

2

The most likely answer here is that the UTM initialization is incorrect. Universal Transverse Mercator is a family of projections and by not specifying the correct zone the pyproj assumes central meridian at -183°: Proj('+proj=tmerc +lat_0=0 +lon_0=-183 +k=0.9996 +x_0=500000 +y_0=0 +datum=NAD83 +units=m +no_defs', preserve_units=True) Which is likely why ...

2

The fundamental problem is that the UTM zone is so far off the correct one that the projections don't behave as expected. The utm zone should be 31. It is odd that the central meridian is correct even in this case, but it is. Using the correct zone gives prompt> echo 3 0 | proj +proj=utm +zone=31 +ellips=wgs84 500000.00 0.00 prompt> echo 3 -.0001 |...

2

The issue you are having is due to axis order changes in PROJ 6+: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 One solution is to use always_xy: from shapely.geometry import LineString from shapely.ops import transform from functools import partial import pyproj roadsegment = LineString([(5.318945751388698, 50....

2

The +init= syntax is deprecated. So all you need is the ESRI:102001 part. See: https://pyproj4.github.io/pyproj/stable/gotchas.html#init-auth-auth-code-should-be-replaced-with-auth-auth-code inputGDF.crs = 'esri:102001'

1

UPDATE: After thinking about it, the most efficient method for you to transform the coordinates is probably to not use apply but to use the column array. from pyproj import Proj pp = Proj(proj='utm',zone=10,ellps='WGS84', preserve_units=False) xx, yy = pp(My_data["LON"].values, My_data["LAT"].values) My_data["X"] = xx My_data["Y"] = yy Using Transformer ...

1

You are on the right track. However, it isn't exactly the same as the proj parameter. There are some differences. But, for your purposes of checking if it is the geos projection, the coordinate operation should provide what you need. >>> from pyproj import CRS >>> geos_dict = {'a': '6378137.0', 'b': '6356752.31414', 'h': '35786023.0', '...

1

From projinfo application of the PROJ library: c:\>projinfo -s EPSG:4979 -t EPSG:7417 returns: Operation n1: unknown id, Inverse of DVR90 height to ETRS89 + Inverse of ETRS89 to WGS 84 (1) + UTM zone 33N, unknown accuracy, Denmark - onshore, at least one grid missing PROJ string: +proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=...

1

There is two way to solve your problem: 1st you can try to use the EPSG code for the EGM2008 height with pyproj which is EPSG:3855 (instead of the NAD83 datum that you are trying to use because it is a horizontal datum and not a vertical datum) 2nd you use PROJ6 instead of pyproj (pyproj is the python library from proj) You need to use the proj ...

1

The issue you are running into is a pyproj 1x issue. I would recommend installing pyproj 2+. You can do this with conda using the conda-forge channel: conda config --add channels conda-forge conda config --add channels conda-forge conda config --set channel_priority strict conda create -n geo geopandas conda activate geo Or, you can install the latest ...

1

The issue you are running into is a pyproj 1x issue. I would recommend installing pyproj 2+. You can do this with conda using the conda-forge channel: conda config --add channels conda-forge conda config --add channels conda-forge conda config --set channel_priority strict conda create -n geo geopandas conda activate geo Or, you can install the latest ...

1

If you have an accuracy near centimeter-level and you are wondering if you need to convert from WGS84 to GRS80 ellipsoid height, there is no reason to The difference between both ellipsoids is at a sub-millimeter level (US and European made a different rounding on the ellipsoid flattening) WGS84 doesn't have any geoid definition. The use of ellipsoid or ...

1

I believe I have now resolved this. The polygon in the .kml file is missing a closing point (i.e. the same as the first point). Adding the first point at the end of the coordinates tag, importing to Google Earth Pro and measuring the area produces 955km^2, much closer to that produced by the python code.

1

So @Spacedman pointed me in the right direction, all praise goes to him. I just want to show the complete code for my question: import pyproj as proj ... crs_wgs = proj.Proj(init='epsg:4326') # assuming you're using WGS84 geographic #Erect own local flat cartesian coordinate system cust = proj.Proj("+proj=aeqd +lat_0={0} +lon_0={1} +datum=WGS84 +units=m"....

1

You can do this but I'm not sure its a good thing. I'd work in EPSG:27700 because all my other map data is most likely to have this, and if I lose the projection info for something with a custom CRS then guessing it is very very hard. But here's how to do it: Suppose we start at 3W, 58N. What's that in EPSG:27700? >>> proj.transform(crs_wgs, ...

1

According to the website, the first box is calculating the distance with the Great-circle distance between two point. This is the geodesic calculation in PROJ (https://proj.org/geodesic.html). To perform this calculation in pyproj, you can use pyproj.Geod.inv(). The distances are in meters. >>> from pyproj import Geod >>> gg = Geod(ellps="...

1

From the pyproj FAQ page, it recommends using the Transformer class for this type of operation as it handles datum shifts. I would recommend reading http://pyproj4.github.io/pyproj/stable/gotchas.html#proj-not-a-generic-latitude-longitude-to-projection-converter as it will help clarify if you should use the geodetic CRS or EPSG:4326 in your use case. So, ...

1

I had a hard time tracking this down but finally looked in the pyproj tests for CRS. What you need to do is test whether the input CRS is geographic or not. Test is directly from pyproj python script (above link). def test_is_geographic(): assert CRS({"init": "EPSG:4326"}).is_geographic is True assert CRS({"init": "EPSG:3857"}).is_geographic is ...

1

Based on your input bounds, it appears that the latlon bounds are likely a mask. The shape of the latlon bounds are not rectangular in the latlon form: Or in the projected form: Also, if you create a transform from the bounds: The xscale and yscale are ~2.6 and ~1.3 respectively. Which is much smaller than the expected 2k-3k. So, the data is likely a ...

1

Since I still want to use pyproj 1.9.6 because Pyproj 2.1.3. still hiccup with the packages that depend on Qt (matplotlib, Jupyter, etc). https://github.com/pyproj4/pyproj/issues/268 I had the same issue, and came across these pages: github.com/OpenDrift/opendrift/issues/94 github.com/OSGeo/proj.4/issues/1191 Basically, it seems to be an issue in version ...

Only top voted, non community-wiki answers of a minimum length are eligible