[gdal-dev] reproject points using raster info

[Frank Warmerdam]

Yann Chemin wrote:
> Hello,
> 
> (C API)
> 
> The aim is to get a (lat,long) pair
> into (projX,projY) pair from hDataset1 projection info.
> 
> Lat/long projection was defined with:
> 
> char    *pszProj4 = "+proj=ll +ellps=wgs84 +datum=wgs84";

Yann,

For a geographic (long/lat) coordinate system you should be
using +proj=longlat, not +proj=ll.

> OGRSpatialReferenceH hSRS;
> OSRImportFromProj4(hSRS, pszProj4);
> printf( "%s\n", pszProj4 );

The Import call looks fine.

> Using GDALGetProjectionRef( hDataset1 ),
> the raster projection info is extracted.
> 
> Here somehow the extracted projection information
> should be imported into a second OGRSpatialReferenceH,
> called hSRS1 for example. But I don't know how to do that.


> I am confused by the C++ examples actually,
> and have difficulties to find a proper C example somewhere to do that.
> 
> thanks for any help/suggestion,
> Yann
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
> 

You might find the code from gdalinfo.c which does roughly the opposite
(projected dataset corners to lat/long) helpful.  Some key parts are:

/* -------------------------------------------------------------------- */
/*      Setup projected to lat/long transform if appropriate.           */
/* -------------------------------------------------------------------- */
     if( GDALGetGeoTransform( hDataset, adfGeoTransform ) == CE_None )
         pszProjection = GDALGetProjectionRef(hDataset);

     if( pszProjection != NULL && strlen(pszProjection) > 0 )
     {
         OGRSpatialReferenceH hProj, hLatLong = NULL;

         hProj = OSRNewSpatialReference( pszProjection );
         if( hProj != NULL )
             hLatLong = OSRCloneGeogCS( hProj );

         if( hLatLong != NULL )
         {
             CPLPushErrorHandler( CPLQuietErrorHandler );
             hTransform = OCTNewCoordinateTransformation( hProj, hLatLong );
             CPLPopErrorHandler();

             OSRDestroySpatialReference( hLatLong );
         }

         if( hProj != NULL )
             OSRDestroySpatialReference( hProj );
     }

...
/************************************************************************/
/*                        GDALInfoReportCorner()                        */
/************************************************************************/

static int
GDALInfoReportCorner( GDALDatasetH hDataset,
                       OGRCoordinateTransformationH hTransform,
                       const char * corner_name,
                       double x, double y )

{
     double	dfGeoX, dfGeoY;
     double	adfGeoTransform[6];

     printf( "%-11s ", corner_name );

/* -------------------------------------------------------------------- */
/*      Transform the point into georeferenced coordinates.             */
/* -------------------------------------------------------------------- */
     if( GDALGetGeoTransform( hDataset, adfGeoTransform ) == CE_None )
     {
         dfGeoX = adfGeoTransform[0] + adfGeoTransform[1] * x
             + adfGeoTransform[2] * y;
         dfGeoY = adfGeoTransform[3] + adfGeoTransform[4] * x
             + adfGeoTransform[5] * y;
     }

     else
     {
         printf( "(%7.1f,%7.1f)\n", x, y );
         return FALSE;
     }

/* -------------------------------------------------------------------- */
/*      Report the georeferenced coordinates.                           */
/* -------------------------------------------------------------------- */
     if( ABS(dfGeoX) < 181 && ABS(dfGeoY) < 91 )
     {
         printf( "(%12.7f,%12.7f) ", dfGeoX, dfGeoY );

     }
     else
     {
         printf( "(%12.3f,%12.3f) ", dfGeoX, dfGeoY );
     }

/* -------------------------------------------------------------------- */
/*      Transform to latlong and report.                                */
/* -------------------------------------------------------------------- */
     if( hTransform != NULL
         && OCTTransform(hTransform,1,&dfGeoX,&dfGeoY,NULL) )
     {

         printf( "(%s,", GDALDecToDMS( dfGeoX, "Long", 2 ) );
         printf( "%s)", GDALDecToDMS( dfGeoY, "Lat", 2 ) );
     }

     printf( "\n" );

     return TRUE;
}


-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | Geospatial Programmer for Rent