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["EPSG","6261"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["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 ); it returns both the lon and lat as -infinity (i.e., the double representation of negative infinity). I don't understand why it's returning invalid values like that. I believe I've used pj_transform() with other lat/lon projections (going from these other lat/lon projections to WGS84 as above) and it's acted basically as a no-op. However, my concern is that it's not quite a no-op, but may slightly change the lat/lon values to account for differences in the datum and so on. So, I *could* just not do the pj_transform() when pj_is_latlong( sourceProjection ) == 1, but it seems like I *do* want to do transform in general, because of datum differences and so on. So my question is: (a) Is it misguided to try to call pj_transform() to go from one lat/lon projection (as in this "Merchich" projection) to a WGS84? (b) If not, do you know why it seems to be failing? I found this: http://trac.osgeo.org/osgeo4w/ticket/320. But I don't believe it's related. Thanks for any help.
_______________________________________________ gdal-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/gdal-dev
