Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I boiled the problem down to that python example:

import sys
try:
from osgeo import gdal
from osgeo import osr
except:
import gdal
import osr

if __name__ == '__main__':
version_num = int(gdal.VersionInfo('VERSION_NUM'))
print( "GDAL Version: " + str(version_num))

wkt = """
PROJCS["Popular Visualisation CRS / Mercator",
    GEOGCS["Popular Visualisation CRS",
        DATUM["Popular_Visualisation_Datum",
            SPHEROID["Popular Visualisation Sphere",6378137,0,
                AUTHORITY["EPSG","7059"]],
                TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6055"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4055"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
    AUTHORITY["EPSG","9001"]],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],
    AUTHORITY["EPSG","3785"]]"""

X = 2662981.191207001451403
Y = -1935374.062384711578488

srs = osr.SpatialReference()
srs.ImportFromWkt(wkt)

srsLatLong = srs.CloneGeogCS()
ct = osr.CoordinateTransformation(srs, srsLatLong)
(int, lat, height) = ct.TransformPoint(X, Y)

# Report results
print('X: %s\t\tY: %s' % (X, Y))
print('Lat: %s\t\tLon: %s' % (gdal.DecToDMS(lat, 'Lat', 2), gdal.DecToDMS(int, 'Long', 2)))

my output is:

GDAL Version: 1900
X: 2662981.19121                Y: -1935374.06238
Lat:  17d 1' 1.01"S             Lon:  23d55'19.08"E

However, when I run this from proj4, I get a different answer:

invproj +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs
2662981.191207001451403 -1935374.062384711578488
23d55'19.081"E  17d7'29.797"S

The latitude is off by 6'28" or so... ?????? What could be going wrong? Many thanks, Pat.

share|improve this question
    
Here is the problem I found. When you use gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3785... the source projection is changed from 4326 to 7059... this makes the invproj not as expected. Use newer EPSG:3857 instead. –  Pat Cappelaere Nov 29 '12 at 14:28

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.