3

I have following problem:

I want to convert EPSG 3035 coordinates to GPS Latitude & Longitude coordinates (EPSG 4326) via Python.

Therefore I'm using GDAL and following Python Code with EPSG 3035 coordinates: N2628/E4704

pointX = 2628 pointY = 4704 # Spatial Reference System inputEPSG = 3035 outputEPSG = 4326 # create a geometry from coordinates point = ogr.Geometry(ogr.wkbPoint) point.AddPoint(pointX, pointY) print pointX, pointY # create coordinate transformation inSpatialRef = osr.SpatialReference() inSpatialRef.ImportFromEPSG(inputEPSG) outSpatialRef = osr.SpatialReference() outSpatialRef.ImportFromEPSG(outputEPSG) coordTransform = osr.CoordinateTransformation(inSpatialRef, outSpatialRef) # transform point point.Transform(coordTransform) # print point in EPSG 4326 print point.GetX(), point.GetY()

unfortunately this returns me a a point somewhere in the South Atlantic, but it should be somewhere in Austria.

The Data Vendor of the EPSG 3035 Coordinates mentions following Reference System of the layer:

+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

What am I doing wrong?

Improve this question
9
  • I know this might sound condescending, but rest assure I made this mistake often enough: have you checked that you didn’t mix up X and Y? Plus: could you please give us an example coordinate pair? EDIT: sorry, just saw there is an example coordinate pair Commented Jan 31, 2017 at 17:30
  • Hi Christoph, I've also tried it with swapping the coordinates -> maybe an interesting hint, that it results in the same point in South Atlantic Commented Jan 31, 2017 at 17:37
  • I tried that, too, now … What I find odd, is that your input coordinates are outside of the bounds of the projection (compare spatialreference.org/ref/epsg/3035 ) – I think that’s the actual culprit Commented Jan 31, 2017 at 17:39
  • Maybe, but these are exactly the Coordinates I get from my Data Vendor, I've also tried to add between one and three 0s at the end, unfortunately this doesn't bring me to Austria either :( Commented Jan 31, 2017 at 17:45
  • if you add x_0 and y_0 (false easting, false northing), you at least end up in Europe: coordTransform.TransformPoint(4321000+4704,3210000+2628) >>> (10.068529554343733, 52.02359886909443, 0.0) Commented Jan 31, 2017 at 17:49

2 Answers 2

Reset to default
1

The center of Vienna, in EPSG:3035 coordinate system, is 4,794,100 x 2,808,600 So it seems you have both X and Y reversed, and missing three digits.

Using the proj4 cs2cs utility with your coords in the correct order I get:

echo "4704000 2628000" | cs2cs +init=epsg:3035 +to +init=epsg:4326 15d0'8.974"E 46d38'57.095"N 0.000

(south west of Graz)

Improve this answer
0

Multiply by 1000 and get the x-y order correct and you end up near Vienna. Here's it done in R:

Make a data frame with the coordinates and an id:

> pt2 = data.frame(x=4704000,y=2628000,i=1)

Make it spatial by telling it what its coordinates and projection are:

> coordinates(pt2)=~x+y > proj4string(pt2)="+init=epsg:3035"

Now transform to epsg 4326 long-lat:

> spTransform(pt2,"+init=epsg:4326") coordinates i 1 (15.00249, 46.64919) 1

Here's Vienna in long-lat:

> coordinates(vienna) lng lat 1 16.37208 48.20849

Close enough to still be in Austria? I can test that with my geonames interface functions:

> GNcountryCode(lng=15.00249, lat=46.64919) $languages [1] "de-AT,hr,hu,sl" $distance [1] "0" $countryCode [1] "AT" $countryName [1] "Republic of Austria"

I get the same result using your python code:

>>> pointX=4704000 >>> pointY=2628000 [...] >>> print point.GetX(), point.GetY() 15.0024926758 46.6491931687
Improve this answer
1
  • and that's how it's done ! THX a lot Spacedman, I've tried to add 3 zeros before, and also swapped the order, but I've never tried both at the same time. Commented Jan 31, 2017 at 18:18

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.