Another pixel to lat,lon question

Question

Problem

I have a screenshot of a section of a map (like a tile) in a numpy array.

• The projection of the map is UTM.
• I know the width and height in pixels of the image.
• I do not know the total pixel width and height of the assembled tiles for the whole world.
• I also know the corresponding lat/longs for each corner of the tile, as well as 10 (I can get more if necessary) GCPs throughout the image.
• The image is always north up with no rotation.

I need a python function that accepts my gcp list, the image, and an arbitrary pixel x,y and returns the longitude and latitude for the x,y pixel pair. I have found many related questions but none seem to get me there. I am trying to use gdal but am open to any modern python library since I have to do this many times.

Attempted Solution

Here is what I have so far:

1. Make numpy array gdal friendly:

def array2Gtiff(sct, array):
    cols = array.shape[1]
    rows = array.shape[0]
    outRaster = driver.Create('test.tiff', cols, rows, 1, gdal.GDT_Byte)
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array[:,:,0])
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(3395)
    references = getReferences(sct)
    gcps = []
    for r in references:
        lat = r[0]
        y = r[1]
        lon = r[2]
        x = r[3]
        gcps.append(gdal.GCP(float(lon), float(lat), 0, float(x), float(y)))

    outRaster.SetGCPs(gcps, outRasterSRS.ExportToWkt())
    error_threshold = 0.125
    resampling=gdal.GRA_NearestNeighbour
    tmp_ds = gdal.AutoCreateWarpedVRT( outRaster,
                                  None, 
                                   outRasterSRS.ExportToWkt(),
                                   resampling,
                                   error_threshold)
    return (outputRaster, tmp_ds)

2. Try to get long and lat:

def getPositionOfSymbolUsingGDAL(rasters, symbol):
    old = rasters[0].GetProjectionRef()
    new = rasters[1].GetProjectionRef()

    x = symbol[0]
    y = symbol[1]

    t = osr.CoordinateTransformation(old, new)
    lat, lon, z = t.TransformPoint(x, y)
    return (lon, lat)
    
 

Edit

Okay attempt 2

def array2Gtiff(sct, array):
    
    cols = array.shape[1]
    rows = array.shape[0]
    driver = gdal.GetDriverByName('MEM')
    outRaster = driver.Create('test.tiff', cols, rows, 1, gdal.GDT_Byte)
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array[:,:,0])
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(4326)

    
    references = getReferences(sct)
    gcps = []
    for r in references:
        lat = r[0]
        y = r[1]
        lon = r[2]
        x = r[3]
        gcps.append(gdal.GCP(float(lon), float(lat), 0, float(x), float(y)))

    outRaster.SetGCPs(gcps,  outRasterSRS.ExportToWkt())
  
    error_threshold = 0.01
    resampling=gdal.GRA_NearestNeighbour
    tmp_ds = gdal.AutoCreateWarpedVRT( outRaster,
                                  None, 
                                   outRasterSRS.ExportToWkt(),
                                   resampling,
                                   error_threshold)
  
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()
    return (outRaster,tmp_ds)

def getPositionOfSymbolUsingGDAL(rasters, symbol):

    
    old = rasters[0].GetProjectionRef()
    new = rasters[1].GetProjectionRef()

    xoffset, px_w, rot1, yoffset, rot2, px_h = rasters[1].GetGeoTransform()



    x = symbol[0]
    y = symbol[1]

    posX = px_w * x + rot1 * y + xoffset
    posY = rot2 * x + px_h * y + yoffset

    old_cs= osr.SpatialReference()
    old_cs.ImportFromWkt(old)


    new_cs= osr.SpatialReference()
    new_cs.ImportFromWkt(new)

    t = osr.CoordinateTransformation(old_cs, new_cs)
    lat, lon, z = t.TransformPoint(posX, posY)
    return (lon, lat)

Ok getting super close. This gives lons and lats in the right ballpark but off by 2 or 3 degrees around lat 50 and lon -16. Any other pointers?

Edit 2:

I am closer than ever. Field of view is roughly 7 degrees of longitude (at latitude 52), roughly 1 degree of latitude.

This is the expected list of lons and lats from another but much slower method that I have of obtaining the lon and lat from xy: TRUE LIST: (52.828163, -21.983257) (53.121275, -23.904298) (53.307069, -22.928915) (53.441183, -21.808616) (53.488742, -25.567096) (53.699348, -24.666027) (53.735629, -23.664633) (53.782858, -22.170902) (53.872774, -25.522507)

CALCULATED WITH GDAL: (52.825434609375, -22.060532004780878) (53.115689671874996, -23.943306352988046) (53.300699578125, -22.98735225936255) (53.434748671875, -21.889370700398405) (53.482385953125, -25.57298047450199) (53.693984109375, -24.689860978486056) (53.730542953124996, -23.70841477569721) (53.778180234375, -22.24443936374502) (53.869023421875, -25.529279715936255)

Updated Code:

def array2Gtiff(sct, array):
    # tl and br in (lat, lon)
    tl = getTopLeftPostition(sct)
    br = getBottomRightPosition(sct)
    deltaX =  abs(tl[1] - br[1])
    deltaY = abs(tl[0] - br[0])

    
    cols = array.shape[1]
    rows = array.shape[0]
    
    driver = gdal.GetDriverByName('MEM')
    old = driver.Create('test.tiff', cols, rows, 1, gdal.GDT_Byte)

   
    old.SetGeoTransform((tl[1], deltaX/cols, 0, tl[0], 0, -deltaY/rows))
    outband = old.GetRasterBand(1)
    outband.WriteArray(array[:,:,0])
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(4326)

   
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    
    error_threshold = 0.00001
    resampling=gdal.GRA_NearestNeighbour

    new = gdal.AutoCreateWarpedVRT( old,
                                   None, 
                                   outRasterSRS.ExportToWkt(),
                                   resampling,
                                   error_threshold)
    
    outband.FlushCache()
    return (old, new)

def getPositionOfSymbolUsingGDAL(rasters, symbol):

    # rasters is the output from array2Gtiff
    old = rasters[0].GetProjectionRef()
    new = rasters[1].GetProjectionRef()

    x = symbol[0] #The xth pixel of the image - I want the lon equivalent of this pixel value.
    y = symbol[1] # The yth pixel of the image - I want the lat equivalent of this pixel value.

    xoff, a, b, yoff, d, e = rasters[0].GetGeoTransform()

    xp = (a * x + b * y + xoff) 
    yp = (d * x + e * y + yoff) 

  
 
    old_cs= osr.SpatialReference()
    old_cs.ImportFromWkt(old)


    new_cs= osr.SpatialReference()
    new_cs.ImportFromWkt(new)

    t = osr.CoordinateTransformation(old_cs, new_cs)
    lat, lon, z = t.TransformPoint(xp, yp)
    return (lon, lat)

This is where I am confused. Right now I feel like I'm using the pixels as the source coordinates. Is this the same thing as just linear interpolation? I no longer use GCP and just use the geotransform.

Will I get more accurate values if I convert to meters and use UTM? Why would using meters instead of pixels change the accuracy of the conversion?

Answer 1

Your image may not be in UTM. Your procedure isn't making any use of that information--it is assigning/forcing the lat/lon EPSG:4326 coordinate system to your image and setting the relative location and sizes of the pixels based on the upper left and lower right corners. A problem this approach could cause is if the pic is rotated or stretched one way or the other, like a UTM projected image would be relative to EPSG:4326, the dx,dy would be stretched/compressed and perhaps nonlinear.

What you have now essentially sets you up to interpolate lat/lon between the UL/LR pixels, (while also packaging some unused GCPs in the file). You could use something like gdalwarp https://gdal.org/programs/gdalwarp.html to make use the GCPs, ignore the UTM information and warp the image to match the EPSG:4326 system with something like gdalwarp -t_srs EPSG:4326 input.tif output.tif

If you know the source is in a particular UTM zone and that the pixel grid is aligned with the UTM coordinate axes, you can use something like the assumptions in the first routine, but you'd need to convert the corner points from lat/lon to the UTM coordinates, then add an intermediate step to warp the (now geolocated) UTM image into a (geolocated) lat/lon image.

The CRSs specify a coordinate system, origin, units, and the mechanics of the coordinate system so that one can translate between systems. Your initial "coordinates" are your image's row and column indices, which if you know it was an image of a UTM projection, are almost the scaling and translation you are calculating. The difference is that you are currently specifying the units as lat/lon degrees from the prime meridian/equator origin rather than in meters from the particular UTM zone's origin.

One trick you might be able to use is that since you know the initial image is UTM based, (a meter-based cylindrical projection wrapped around a line of meridian) you could define your own Transverse Mercator coordinate system to match a corner of your image with a proj4 string like:

outRasterSRS.ImportFromProj4(f"+proj=tmerc +lat_0={OriginY} +lon_0={OriginX}")

...but you'd have to calculate your pixel sizes in meters and change your origin for:

outRaster.SetGeoTransform((0, pixelWidth_m, 0, 0, 0, pixelHeight_m))

Further on your additional question:

Your initial CRS helps translate the origin, units and rotation of your data into a specific geographic position. If you want to measure centimeters from a reference point on an architectural plan, you can define a CRS with the geographical lat/lon location of the reference point as the origin, rotated to align with the plan, and measuring in units of centimeters. If you can locate your data in meters east and meters north from a UTM origin, you can use those coordinates. If you have geographic coordinates, there might be some warping to account for. If you know what coordinate plane represents your data, you can choose a CRS to represent that system (origin, units, and projection)