Problem converting EPSG 3035 to EPSG 4326 with GDAL & Python

Question

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?

Answer 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)

Answer 2

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