[gdal-dev] Converting coordinates (easting/Northing to WGS84 Lat/Long)

Tue Aug 4 09:33:35 EDT 2009

I am still having problems using the geo-transformation utilities in
GDAL.

I don't suppose anyone managed to identify why the following code always
refuses to create a OGRCoordinateTransformation for me? 

The code below is an excerpt from code I use to open geotiffs and
extract geocoords.  However, although I can extract geocoords as
Eastings/Northings, I cannot then convert them to WGS84 lat/long
(essential for my purpose).  I am running this on linux SUSE-10 SP2.

Can anyone help?

// g++ -Wno-deprecated -m64 -g -o georef georef.c++ -lgdal -lm

#include "/usr/local/src/FWTools-2.0.6/include/gdal.h"
#include "/usr/local/src/FWTools-2.0.6/include/gdal_priv.h"
#include "/usr/local/src/FWTools-2.0.6/include/ogr_spatialref.h"
/* Following output is from the MapServer tutorial
...data/raster/shdrlfi020l_ugl.tif file.  This file opens
   fine in 'openev' and gives easting/northing readouts OK.
Driver: GTiff/GeoTIFF
Projection is `'
Origin = (176371.000000,548761.000000)
Pixel Size = (1000.000000,-1000.000000)
ERROR 1: No PROJ.4 translation for source SRS, coordinate
transformation initialization has failed.
!!!Cannot create transform from input file georef space to WGS84
*/
/*  Following output is from a geotiff file that also works fine with
'openev'.
Driver: GTiff/GeoTIFF
Projection is `PROJCS["TM
D-01",GEOGCS["NAD27",DATUM["North_American_Datum_1927",SPHEROID["Clarke
1866",6378206.4,294.9786982139006,AUTHORITY["EPSG","7008"]],AUTHORITY["E
PSG","6267"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AU
THORITY["EPSG","4267"]],PROJECTION["Transverse_Mercator"],PARAMETER["lat
itude_of_origin",0],PARAMETER["central_meridian",-124.5],PARAMETER["scal
e_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],
UNIT["metre",1,AUTHORITY["EPSG","9001"]]]'
Origin = (138586.280052,3768995.081326)
Pixel Size = (50.000000,-50.000000)
ERROR 1: No PROJ.4 translation for source SRS, coordinate
transformation initialization has failed.
!!!Cannot create transform from input file georef space to WGS84
*/
int main(int argc, char **argv)
{
  GDALDataset  *poDataset;

  GDALAllRegister();

  poDataset = (GDALDataset *) GDALOpen( argv[1], GA_ReadOnly );
  if( poDataset == NULL )
  {
    fprintf(stderr, "!!!poDataset was NULL\n");
    exit(-1);
  }
  double adfGeoTransform[6];

  printf( "Driver: %s/%s\n",
          poDataset->GetDriver()->GetDescription(), 
          poDataset->GetDriver()->GetMetadataItem( GDAL_DMD_LONGNAME )
);

  if( poDataset->GetProjectionRef()  != NULL )
      printf( "Projection is `%s'\n", poDataset->GetProjectionRef() );

  if( poDataset->GetGeoTransform( adfGeoTransform ) == CE_None )
  {
      printf( "Origin = (%.6f,%.6f)\n",
              adfGeoTransform[0], adfGeoTransform[3] );

      printf( "Pixel Size = (%.6f,%.6f)\n",
              adfGeoTransform[1], adfGeoTransform[5] );
  }

  // We have to convert from native input file 
  // geospace to our standard WGS84 space so create a transform here.

  // First create a spatial reference in the input space
  OGRSpatialReference iref;
  const char *proj_ref = poDataset->GetProjectionRef();
  char *tmp = NULL;
  if (proj_ref)
  {
    tmp = (char *)malloc(strlen(proj_ref)+1);
    iref.importFromWkt(&tmp);

    // Create a spatial reference in WGS88 lat/long space and
    // create a transform object betseen input and output space
    OGRSpatialReference oref;
    oref.SetWellKnownGeogCS("WGS84");
    OGRCoordinateTransformation *coord_transform
      = OGRCreateCoordinateTransformation(&iref, &oref);
    if (!coord_transform)
    {
      fprintf(stderr, "!!!FAILED! Cannot create transform from input
file "
              "georef space to WGS84\n");

    } // end case default
    free(tmp);

  } // end if got projection ref

} // end main

Gary 

________________________________

From: Smart, Gary 
Sent: 21 July 2009 08:46
To: 'gdal-dev at lists.osgeo.org'
Cc: 'Martin Chapman'
Subject: RE: [gdal-dev] Converting raster pixel space to geospace

Yes I have the gdalinfo - it seems happy enough with the files I am
experimenting with.  I'm not sure what 'proj' is supposed to do and
can't find any info on it - but it is in my path and runs?

I have attached a small extract from my ever failing code - it compiles
and runs.  Perhaps someone could tell me what's wrong.  It seems to fail
to set up a coord transform for any file, even though openev happily
deals with them.  I can't help feeling that there is something decidedly
dodgy about the way I prime the transform with the 'importFromWkt call.
Perhaps the transform is never getting the transform string at all?

I also noticed that one of the map tif files from the mapserver tutorial
doesn't return any projection info fro mthe getProjectionRef call - will
this cause problems when trying to set up a transformation in future.

The files I used, and the response they give, are included in the
attachment as comments.

Gary 

________________________________

From: Martin Chapman [mailto:chapmanm at pixia.com] 
Sent: 20 July 2009 22:46
To: Smart, Gary
Subject: RE: [gdal-dev] Converting raster pixel space to geospace

Gary,

I tried the projection string in some of my code and it worked fine.  I
stepped through the gdal code and proj4 loaded the transform object ok
too which means the projection string is fine and that your environment
suffers from one of the following:

1.    Proj.dll cannot be found in the path on your machine or is not in
the same directory as gdal.dll and is failing when gdal calls
LoadLibrary() on "proj.dll".

2.   The proj.dll can be found but a dependency like the nad shift files
or gdal data cannot be found.

I would try the following:

1.   Use gdalinfo.exe that should come with fwtools or can be downloaded
from gdal.org and run the command gdalinfo.exe - formats and or
ogrinfo.exe -formats.

2.   See if the output from that command gives any useful information on
why proj is failing.

3.   From the command line type in proj.exe and see if the proj util can
be found in your path.  If it can be found then maybe it will give some
useful info.

Sorry I am not more help than that but the projection string is valid.

Martin

From: gdal-dev-bounces at lists.osgeo.org
[mailto:gdal-dev-bounces at lists.osgeo.org] On Behalf Of Smart, Gary
Sent: Monday, July 20, 2009 8:50 AM
To: gdal-dev at lists.osgeo.org
Subject: RE: [gdal-dev] Converting raster pixel space to geospace

Thanks Martin, but now to the next problem...

I downloaded PRROJ.4 which seemed to configure and compile OK.

Having set the PROJ_LIB as you proposed, I now get:

ERROR 1: No PROJ.4 translation for source SRS, coordinate transformation
initialization has failed.

I guess I am still missing comething.  The failure came from the same
geotiff file as below.

I also couldn't find a 'data dir' in the GDAL source and therefore
GDAL_DATA env var remained unset.  We installed everything from the
FWTools download for our Linux SUSE10 system.  Should there be a 'data'
directory under the FWTools directory somewhere?  There is no GDAL
source in this directory - just include files and libraries.  I presumed
this was all that was needed to use GDAL?   

________________________________

From: Martin Chapman [mailto:chapmanm at pixia.com] 
Sent: 20 July 2009 14:49
To: Smart, Gary; gdal-dev at lists.osgeo.org
Subject: RE: [gdal-dev] Converting raster pixel space to geospace

Gary,

Set your paths like such on startup or as an environment variable:

GDAL_DATA points at the data dir that comes with gdal source.

PROJ_LIB point to the nad directory that you download from proj4
website.

string env = string("GDAL_DATA=") + dataPath;

_putenv(env.c_str());

env = string("PROJ_LIB=") + nadPath;

_putenv(env.c_str());

and / or

CPLSetConfigOption("GDAL_DATA", dataPath.c_str());

CPLSetConfigOption("PROJ_LIB", nadPath.c_str());

Martin

From: gdal-dev-bounces at lists.osgeo.org
[mailto:gdal-dev-bounces at lists.osgeo.org] On Behalf Of Smart, Gary
Sent: Monday, July 20, 2009 7:38 AM
To: gdal-dev at lists.osgeo.org
Subject: RE: [gdal-dev] Converting raster pixel space to geospace

OK - I can get it to compile and not crash instantly (bonanza!)

However, when I try to create the OGR transform using importFromWkt, 

I get the following:

ERROR 6: Unable to load PROJ.4 library (libproj.so), creation of

OGRCoordinateTransformation failed.

Does anyone know what this means?  Is there some part of OGR/GDAL I have
yet to install?

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20090804/71041e6b/attachment-0001.html