GIS / GEOTiff / GDAL / Python How to get coordinates from pixel

Question

I'm working on the project to detect the object from GEOTiff files and return coordinates of the objects and those output will use for drone to fly to those coordinate

I use tensorflow with YOLO v2(image detector framework) and OpenCV to detect the objects that I need in GEOTiff

import cv2
from darkflow.net.build import TFNet
import math
import gdal

# initial stage for YOLO v2 
options = {
    'model': 'cfg/yolo.cfg',
    'load': 'bin/yolov2.weights',
    'threshold': 0.4,
}
tfnet = TFNet(options)

# OpenCV read Image
img = cv2.imread('final.tif', cv2.IMREAD_COLOR)
img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)

#Predict the image
result = tfnet.return_predict(img)

#Calculate the center and radius of each objects
i = 0
while i < len(result):
    tl = (result[i]['topleft']['x'], result[i]['topleft']['y'])
    br = (result[i]['bottomright']['x'], result[i]['bottomright']['y'])
    point = (int((result[i]['topleft']['x']+result[i]['bottomright']['x'])/2), int((result[i]['topleft']['y']+result[i]['bottomright']['y'])/2))
    radius = int(math.hypot(result[i]['topleft']['x'] - point[0], result[i]['topleft']['y'] - point[1]))
    label = result[i]['label']
    result[i]['pointx'] = point[0]
    result[i]['pointy'] = point[1]
    result[i]['radius'] = radius    
    i += 1

print(result)

So the results come out like set of JSON

[{'label': 'person', 'confidence': 0.6090355, 'topleft': {'x': 3711, 'y': 1310}, 'bottomright': {'x': 3981, 'y': 1719}, 'pointx': 3846, 'pointy': 1514, 'radius': 244}]

as you can see the location of the object is return in pixel (x,y) and I want to use these x,y to convert to coordinate in lat,lng so I try to use GDAL (the library use for read the GEO infomation that contain inside the image)

so here's the GEO infomation of the image by using gdalinfo in terminal

Driver: GTiff/GeoTIFF
Files: final.tif
Size is 8916, 6888
Coordinate System is:
PROJCS["WGS 84 / UTM zone 47N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",99],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32647"]]
Origin = (667759.259870000067167,1546341.352920000208542)
Pixel Size = (0.032920000000000,-0.032920000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_SOFTWARE=pix4dmapper
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  667759.260, 1546341.353) (100d33'11.42"E, 13d58'57.03"N)
Lower Left  (  667759.260, 1546114.600) (100d33'11.37"E, 13d58'49.65"N)
Upper Right (  668052.775, 1546341.353) (100d33'21.20"E, 13d58'56.97"N)
Lower Right (  668052.775, 1546114.600) (100d33'21.15"E, 13d58'49.59"N)
Center      (  667906.017, 1546227.976) (100d33'16.29"E, 13d58'53.31"N)
Band 1 Block=8916x1 Type=Byte, ColorInterp=Red
  NoData Value=-10000
Band 2 Block=8916x1 Type=Byte, ColorInterp=Green
  NoData Value=-10000
Band 3 Block=8916x1 Type=Byte, ColorInterp=Blue
  NoData Value=-10000
Band 4 Block=8916x1 Type=Byte, ColorInterp=Alpha
  NoData Value=-10000

Any one?

Parker

Answer 1

You need to transform pixel coordinates to geographic space using the GeoTransform matrix that is associated to your raster files. Using GDAL you could do something like the following:

# open the dataset and get the geo transform matrix
ds = gdal.Open('final.tif') 
xoffset, px_w, rot1, yoffset, px_h, rot2 = ds.GetGeoTransform()

# supposing x and y are your pixel coordinate this 
# is how to get the coordinate in space.
posX = px_w * x + rot1 * y + xoffset
posY = rot2 * x + px_h * y + yoffset

# shift to the center of the pixel
posX += px_w / 2.0
posY += px_h / 2.0

Of course the position you get will be relative to the same coordinate reference system that is used for your raster dataset. So if you need to transform it to lat/long you will have to do further elaborations:

# get CRS from dataset 
crs = osr.SpatialReference()
crs.ImportFromWkt(ds.GetProjectionRef())
# create lat/long crs with WGS84 datum
crsGeo = osr.SpatialReference()
crsGeo.ImportFromEPSG(4326) # 4326 is the EPSG id of lat/long crs 
t = osr.CoordinateTransformation(crs, crsGeo)
(lat, long, z) = t.TransformPoint(posX, posY)

Sorry I'm not really fluent in python, so probably you will have to adapt this code. Checkout the documentation of GeoTransform here for the C++ API to learn more about the matrix elements.

Novercal

Great example and followup. It appears the return values from GetGeoTransform have changed to xoffset, px_w, rot1, yoffset, rot2, px_h from xoffset, px_w, rot1, yoffset, px_h, rot2 – Toplevel

Thanks, this was helpful to me along with user2897775's comment. If I'm not mistaken, the order of return values in the last line should be (long, lat, z) = t.TransformPoint(posX, posY). – Argile

Answer 2

Without the excellent and clear Python code posted by Gabriella, I don't know if I would have ever figured out how to do this in C. The documentation and examples for gdal are amazingly sparse.

Here's a C version of Gabriella's code:

const char fn[] = "/path/to/geo/file.tif";

GDALDatasetH  hDataset;
GDALAllRegister(); // Register all GDAL formats
hDataset = GDALOpen( fn, GA_ReadOnly ); // Open our geo file (GeoTIFF or other supported format)
if (hDataset == NULL)
{
    printf("Failed to open dataset\n");
    return;
}

// These are the input points to be transformed, in pixel coordinates of the source raster file
double x = 20;
double y = 20;

double        adfGeoTransform[6];
GDALGetGeoTransform( hDataset, adfGeoTransform );

// Put the returned transform values into named vars for readability
double xoffset = adfGeoTransform[0];
double px_w = adfGeoTransform[1];
double rot1 = adfGeoTransform[2];
double yoffset = adfGeoTransform[3];
double rot2 = adfGeoTransform[4];
double px_h = adfGeoTransform[5];

// Apply transform to x,y. Put into posX,posY
double posX = px_w * x + rot1 * y + xoffset;
double posY = rot2 * x + px_h * y + yoffset;

// Transform to center of pixel
posX += px_w / 2.0;
posY += px_h / 2.0;

OGRErr err = 0;

// sr0 is the "from" spatial reference, pulled out of our file
OGRSpatialReferenceH sr0 = OSRNewSpatialReference(GDALGetProjectionRef(hDataset));

// sr1 is the "to" spatial reference, initialized as EPSG 4326 (lat/lon)
OGRSpatialReferenceH sr1 = OSRNewSpatialReference(NULL);
err = OSRImportFromEPSG(sr1, 4326);

double xtrans = posX;
double ytrans = posY;
double ztrans = 0;
int pabSuccess = 0;

// Make our transformation object
OGRCoordinateTransformationH trans = OCTNewCoordinateTransformation(sr0, sr1);

// Transform our point posX,posY, put it into xTrans,yTrans
OCTTransformEx(trans, 1, &xtrans, &ytrans, &ztrans, &pabSuccess);
    
GDALClose(hDataset);

printf("map coordinates (%f, %f)\n", xtrans, ytrans);

Barre