If You Have A GeoTiff, Would It Be Possible To Transform A Lat/Lon Point To An X,Y Using The GeoTransform?

Question

I'm using the GDAL library. Currently, I can take in an upper left point and an upper right point and chip an image out of the original. What I'd like to do now, is take in two WKT points and convert to X,Y coordinatees to do the same thing. I was just wondering if it was possible to do this if I knew the GeoTransform and what coordinate system it was using (WGS84)?

Answer 1

I ran into this before too, and here's a nice enough way of doing coordinate transforms.

Note from GDAL documentation:

The coordinate system returned by GDALDataset::GetProjectionRef() describes the georeferenced coordinates implied by the affine georeferencing transform returned by GDALDataset::GetGeoTransform().

We can use this with OGRCoordinateTransformation to do the transformation for us.

Basically the code will look something like this:

// Load up some dataset.
dataset = (GDALDataset *) GDALOpen( mapfile, GA_ReadOnly );

// Define Geographic coordinate system - set it to WGS84.
OGRSpatialReference *poSRS_Geog = new OGRSpatialReference();
poSRS_Geog->importFromEPSG( 4326 ); // WGS84

// Define Projected coordinate system - set to the GeoTransform.
const char *sProj = dataset->GetProjectionRef();
OGRSpatialReference *poSRS_Proj = new OGRSpatialReference( sProj );

// Set up the coordinate transform (geographic-to-projected).
OGRCoordinateTransformation *poCT_Geog2Proj;
poCT_Geog2Proj = OGRCreateCoordinateTransformation( poSRS_Geog, poSRS_Proj );

// Now everything is set up and we set transforming coordinates!
// Pass Lon/Lat coordinates to the Transform function:
double x = lon;
double y = lat;
poCT_Geog2Proj->Transform( 1, &x, &y );

// Now x and y variables will contain the X/Y pixel coordinates.

That's how you can convert between longitude/latitude and pixel coordinates. Note you can use arrays with Transform(), and convert multiple coordinates together. The first argument is the number of coordinate pairs to transform, and the second and third arguments are pointers to the x's and y's. I just transform one pair here.

Note it's equally easy to set up the inverse transform:

// Set up the coordinate transform (projected-to-geographic).
OGRCoordinateTransformation *poCT_Proj2Geog;
poCT_Proj2Geog = OGRCreateCoordinateTransformation( poSRS_Proj, poSRS_Geog );