0

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)
  1. 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?

6
  • What is this a re-post of? Please provide links to at least some of the "many related questions" that you mention.
    – PolyGeo
    Nov 18, 2020 at 21:35
  • The question looks edited to now assign the EPSG:3395 World Mercator projection (meters east and north from 0N,0W : epsg.io/map#srs=3395&x=0&y=0&z=2&layer=streets ). You need to transform your GCPs from lat,lon (EPSG:4326) to match the assigned projection before you store them as GCPs. Or you leave GCPs in lat/lon, assign EPSG:4326 and depend on the warping algorithm to un-project the screencaptured image into a lat/lon aligned raster.
    – Dave X
    Nov 18, 2020 at 21:57
  • So I assign the initial projection EPSG 4326 and then I'm getting values much much closer to what I expect. Still off by a few degrees, any other pointers? Appreciate your help with this Dave.
    – jweisbaum
    Nov 19, 2020 at 1:32
  • How big are your areas of interest? If the areas are small, the projection warping errors should be small. The best way is to identify the projection/coordinate system of your source data, (UTM? World Mercator?) and then figure out the pixel indices-to-coordinate system's coordinates. From an earlier edit, you already know the lat lons of the corners, so you just need to translate the corner coords from EPSG:4326 to (UTM? ESPG:3395?) and re-apply them with the SetGeoTransform() Per gdal.org/user/raster_data_model.html#affine-geotransform you should do corners or GCPs, not both.
    – Dave X
    Nov 19, 2020 at 3:07
  • For the second getPositionOfSymbolUsingGDAL() function, I'm not sure which are rasters[0] and [1]. If you have the pixel coords in x,y relative to the original, unwarped image, and its (UTM? EPSG:3395?) CRS, do the same translation x,y -> posX,posY in that source CRS (meters from origin), then convert it to EPSG:4326 in the same way to get degrees lon,lat.
    – Dave X
    Nov 19, 2020 at 3:18

1 Answer 1

4

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)

2
  • Thanks so much this is very helpful. I've tried adjusting it to use gdal warp, but the osr.CoordinateTransformation() step fails because my old crs is null. How do I determine what to set it to?
    – jweisbaum
    Nov 18, 2020 at 20:42
  • With distinct files, gdal_translate, gdal.org/programs/gdal_translate.html can assign a CRS. If you know it's UTM in a particular zone like UTM18N, you can set it with an EPSG code like spatialreference.org/ref/epsg/wgs-84-utm-zone-18n . But you also need to set the corner coordinates within that CRS convert your known lat,lons to easting and northing.
    – Dave X
    Nov 18, 2020 at 21:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.