[gdal-dev] geoTransforms and OGRSpatialTransforms

[Smart, Gary]

1)       Can anyone tell me why the following program returns a
GetProjectionRef as 'empty string' (not NULL) from some files and not
others (even though openEv seems to be able to open them and
georeference them all fine?)  I am looking at geotiffs from the
MapServer tutorial directory.

2)       Can anyone offer an explanation as to why - even with a valid
projectionRef, I can never create the OGRSpatialReference

 

I am using Linux SLES10 SP2, gdal-svn-trunk-2009.12.14 and
proj4-4.7.0-1.1 (well - at least I think I am - can someone tell me a
good way to test this?)

 

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

#include <proj_api.h>
#include <gdal.h>
#include <gdal_priv.h>
#include <ogr_spatialref.h>

/* Using from the MapServer tutorial ...data/raster/shdrlfi020l_ugl.tif
file.  
   This file opens fine in 'openev' and gives easting/northing readouts
OK.  
   However - this program fails to manage to create either a pixel->geo
transform OR
   a UTM to Lat/Long transform!  Output from this program is...

Driver: GTiff/GeoTIFF
Input has no ProjectionRef!
No geo transform!
ERROR 1: No PROJ.4 translation for source SRS, coordinate
transformation initialization has failed.
!!!FAILED! Cannot create transform from input space   to WGS84
*/

/* Using the MapServer tutorial
...data/raster/mod09a12003161_ugl_ll_8bit.tif file.
   This file opens fine in 'openev' and gives Lat/Long readouts OK.  
   This program does get a geoTransform but fails to create the 

   UTM to Lat/Long transform too!  Output from this program is...

Driver: GTiff/GeoTIFF
ProjectionRef is 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS
84",6378137,298.2572235630016,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG"
,"6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHOR
ITY["EPSG","4326"]]'
Origin = (-97.374337,49.415919)
Pixel Size = (0.004733,-0.004733)
ERROR 1: No PROJ.4 translation for source SRS, coordinate
transformation initialization has failed.
!!!FAILED! Cannot create transform from input space GEOGCS["WGS
84",DATUM["WGS_1984",SPHEROID["WGS
84",6378137,298.2572235630016,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG"
,"6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHOR
ITY["EPSG","4326"]]  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];

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

  const char *proj_ref = poDataset->GetProjectionRef();
  if( proj_ref != NULL && strlen(proj_ref))
      fprintf(stderr, "ProjectionRef is '%s'\n",
poDataset->GetProjectionRef() );
  else
    fprintf(stderr, "Input has no ProjectionRef!\n");

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

      fprintf(stderr, "Pixel Size = (%.6f,%.6f)\n",
              adfGeoTransform[1], adfGeoTransform[5] );
  }
  else
    fprintf(stderr, "No geo transform!\n");

  // 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;
  char *tmp = NULL;
  if (proj_ref)
  {
    tmp = (char *)malloc(strlen(proj_ref)+1);
    iref.importFromWkt(&tmp);
    //iref.SetWellKnownGeogCS("WGS84");

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

    } // end case default
    free(tmp);

  } // end if got projection ref
} // end main
 

Gary

 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20100107/cef5a12d/attachment.html