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

[Michael Katz - NOAA Affiliate]

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.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20141225/385b78a0/attachment.html>