Getting lat,lon from row, column

28. Getting lat,lon from row, column#

28.1. Introduction#

The function rowcol2latlon defined below takes the name of a tiffile and a row,column pair for the raster and returns the latitude and longitude of that pixel. It uses the pyproj.CRS object to hold the coordinate reference system and the pyproj.Transformer object to do the transform.

from pathlib import Path
from pyproj import Transformer, CRS
import rioxarray

data_dir = Path().home() / 'repos/a301/satdata/landsat'
tif_filename = data_dir / 'vancouver' / f"week4_clipped_TIRS2.tif"
print(f"{tif_filename=}")

tif_filename=PosixPath('/Users/phil/repos/a301/satdata/landsat/vancouver/week4_clipped_TIRS2.tif')

28.2. transforming lat/lon to utm#

Here is how we did this using the code from - Zooming an image.

# zone 10 UTM
proj_code = 32610
#
# we want to go from x,y in utm zone 10 to lat,lon
# first, get x,y from lat, lon, the show you can reverse that
#
p_utm = CRS.from_epsg(proj_code)
p_latlon = CRS.from_epsg(4326)
transform = Transformer.from_crs(p_latlon, p_utm)
ubc_lon = -123.2460 
ubc_lat = 49.2606
ubc_x, ubc_y = transform.transform(ubc_lat, ubc_lon)
print(f"{ubc_x=:.2f} m, {ubc_y:.2f}")

ubc_x=482101.01 m, 5456455.22

#dir(transform)

transform.description

'UTM zone 10N'

28.3. doing the inverse: utm to lat/lon#

To go the other direction, just reverse the parameters in the Transformer call

# zone 10 UTM
proj_code = 32610
#
# we want to go from x,y in utm zone 10 to lat,lon
# first, get x,y from lat, lon, the show you can reverse that
#
p_utm = CRS.from_epsg(proj_code)
p_latlon = CRS.from_epsg(4326)
transform = Transformer.from_crs(p_utm, p_latlon)
ubc_lat, ubc_lon = transform.transform(ubc_x, ubc_y)
print(f"{ubc_lat=:.2f} m, {ubc_lon=:.2f} m")

ubc_lat=49.26 m, ubc_lon=-123.25 m

28.3.1. wrap this in a function#

Find the lat/lon of the upper left corner with row=0, col=0

def rowcol2latlon(tiffile,row,col):
    """
    return the latitude and longitude of a row and column in a geotif

    Parameters
    ----------

    tiffile: Path object
        path to a clipped tiffile
    row: float
       row of the pixel
    col: float
       column of the pixel

    Returns
    -------

    (lon, lat):  (float, float)
       longitude (deg east) and latitude (deg north) on
       a WGS84 datum
    
    """
    has_file = tiffile.exists()
    if not has_file:
        raise IOError(f"can't find {filename}, something went wrong above") 
    the_band = rioxarray.open_rasterio(tif_filename,masked=True)
    epsg_code = the_band.rio.crs.to_epsg()
    p_utm = CRS.from_epsg(epsg_code)
    p_latlon = CRS.from_epsg(4326)
    affine_transform = the_band.rio.transform()
    x, y = affine_transform*(row, col)
    crs_transform = Transformer.from_crs(p_utm, p_latlon)
    lat, lon = transform.transform(x,y) 
    return (lon, lat)

lon, lat = rowcol2latlon(tif_filename,0,0)
print(f"{lon=:.1f} deg,{lat=:.1f} deg")

lon=-123.3 deg,lat=49.3 deg

28.4. Two gotchas to watch out for#

  • Don’t confuse the pyproj transformer.Transform with the affine transform – two completely different things

  • pyproj has an alternate way of doing coordinate transforms using a Proj object and a transform (lower case t) method. These are in lots of examples on the web, but have been superseded by Transformer.transform

Contents