2

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?

  • 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 – christoph Jan 31 '17 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 – Into Numbers Jan 31 '17 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 – christoph Jan 31 '17 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 :( – Into Numbers Jan 31 '17 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) – christoph Jan 31 '17 at 17:49
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
  • 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. – Into Numbers Jan 31 '17 at 18:18
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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.