[gdal-dev] problem projecting "Merchich" to WGS84

[Even Rouault]

Le jeudi 25 décembre 2014 11:02:49, Michael Katz - NOAA Affiliate a écrit :
> I have a certain .tif file that I open with GDALOpen() and then call
> GDALGetProjectionRef(), which returns:
> 
> srs = GEOGCS["Merchich",DATUM["Merchich",SPHEROID["Clarke 1880
> (IGN)",6378249.2,293.4660212936265,AUTHORITY["EPSG","7011"]],AUTHORITY["EPS
> G","6261"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORI
> TY["EPSG","4261"]]
> 
> 
> I make this into a projection with:
> 
> OGRSpatialReference oSRS;
> oSRS.importFromWkt( &srs );
> char *sP = NULL;
> oSRS.exportToProj4( &sP );
> sourceProjection = pj_init_plus( sP );
> CPLFree( sP );
> 
> 
> (here, sP = "+proj=longlat +a=6378249.2 +b=6356515
> +towgs84=31,146,47,0,0,0,0 +no_defs")
> 
> I think make another projection with:
> 
> targetProjection = pj_init_plus( "+proj=latlong +datum=WGS84" );
> 
> 
> I use GDALApplyGeoTransform() to convert the top-left (0, 0) point of the
> .tif file to geo coordinates, which gives me a lat/lon value of (lon = -7,
> lat = 34).
> 
> pj_is_latlong( sourceProjection ) returns 1, which confirms that this file
> has a lat/lon projection, and the (-7, 34) seems to be correct (at least in
> the right area). (The actual numbers from my file are not exact values like
> this, so I'm rounding here, but it shouldn't matter for the purpose of this
> discussion.)
> 
> The problem is that when I call
> 
> pj_transform( sourceProjection, targetProjection, 1, 1, &lon, &lat, NULL );

Michael,

pj_transform expects radians when the input projection is geodetic (and output 
radians when the output projection is geodetic) : 
https://trac.osgeo.org/proj/wiki/ProjAPI

Instead of using directly proj.4, you could use the OSR API that will take 
care of that annoying conversion: http://gdal.org/osr_tutorial.html

Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com