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 used the affine transform for the image to calculate some sample latitudes/longitudes. The only problem I had was if the image is North Facing, geoTransform[2] and geTransform[4] need to be zeroed out when calculating the Lat/Lon.

x = (int)Math.Abs(Math.Round((Latitude - geotransform[0]) / geotransform[1]));
y = (int)Math.Abs(Math.Round((Longitude - geotransform[3]) / geotransform[5]));

If you wanted to brute force it, you could do the following (I did this and it worked but this is just pseudocode):

//Get the Pixel for the length and width, this portion is for the full image
pixelXSize = AbsoluteValue((latitudeAt(Zero)-(latitudeAt(Length)-1))/imageLength);
pixelYSize = AbsoluteValue((longitudeAt(Zero)-(LongitudeAt(Width)-1))/imageWidth);

//Calculate the x,y conversion for the points you want to calculate
x = AbsoluteValue((latitudeToConvert-latitudeAt(Zero))/pixelXSize);
y = AbsoluteValue((longitudeToConvert-longitudteAt(Zero))/pixelYSize);

This answer may be off by one or two pixels. If you are using a geotransform, again the North Facing variables may mess up the answer that is returned. So far, I've only tested this for north facing images.

Answer 2

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 );

Answer 3

I use this method:

void transformCoordinatesEPSG(OGRGeometry &geometry,int from, int to) {
    OGRSpatialReference srcSpatialReference;
    OGRErr error = srcSpatialReference.importFromEPSG(from);

    #ifdef __OGRTRANSFORMDEBUG
        qDebug() << "Import EPSG  " << from << "return " << error;
    #endif

    OGRSpatialReference dstSpatialReference;
    error = error | dstSpatialReference.importFromEPSG(to);

    #ifdef __OGRTRANSFORMDEBUG
    qDebug() << "Import EPSG  " << to << "return " << error;
    #endif

    OGRCoordinateTransformation* coordTrans = OGRCreateCoordinateTransformation(&srcSpatialReference, &dstSpatialReference);
    geometry.transform(coordTrans);
}

For lat/long to must be 4326.