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 _______________________________________________ gdal-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/gdal-dev
