# Can you determine the row,col of an un-orthorectified image corresponding to a lat,lon?

Is there are way to determine the row and column of a pixel in an un-orthorectified image that corresponds to a latitude and longitude given a particular DEM? That is, convert a spatial coordinate to a row,col coordinate for a specific RPC file and DEM source?

Use case: I have a stack of satellite images in nitf format and lat,lon coordinates for points of interest. I want to crop the images around each of the POIs and circumvent any pixel distortions due to orthorectification. So rather than orthorectify the images and crop them around a lat,lon, I want to convert the lat,lon to row,col and crop the un-orthorectified image.

Current approach: The only way I can think to do this for a particular image is to:

1. Create a new 2-band image that is the same size as the original image where the pixel values are the row and column.
2. Orthorectify this new image
3. Sample the new image at a lat,lon to get a row,col coordinate.

For my processing I'm using Python and GDAL bindings.

Problem: This approach seems very inefficient with time and memory, and I don't fully understand the accuracy and precision. So if it's possible, I would prefer to essentially reverse-orthorectify a lat,lon coordinate.

There's a python library for this called orthorectification: https://pypi.org/project/orthorectification/. There is a function called `lon_lat_alt_to_xy` that "Returns an image pixel coordinate (x, y) corresponding to a provided world coordinate (lon, lat, alt) using provided RPC coefficients"

The solution I'm using is to use the forward transform and minimization to find the row, col located at lat, lon of interest.

There is also an analytical form with an example in this package, which is described in this paper section 1.1.

``````import gdal
import numpy as np
from scipy.optimize import minimize

def fun(q, lon, lat, z):
col, row = q
success, pnt = tr.TransformPoint(0, col, row, z)
error = np.sqrt((pnt[0]-lon)**2 + (pnt[1]-lat)**2)
return error

image_path = r'D:\18SEP23102516-M1BS-012454185010_01_P001.NTF'
lat = 59.90823713424692
lon = 19.000549576907073
z = 15 # get by sampling srtm at lat lon

ds = gdal.Open(image_path)
tr = gdal.Transformer(ds, None, ['METHOD=RPC', 'RPC_PIXEL_ERROR_THRESHOLD=0.05'])

# perform minimization
col_0 = 0.5*ds.RasterXSize  # start at center of image
row_0 = 0.5*ds.RasterYSize  # start atcenter of image
res = minimize(fun, (col_0, row_0), args=(lon, lat, z), method='Nelder-Mead')
col, row = res.x
print(col, row)
#  20.499992952535706 10.499958201792301

# go round trip just to check the result against original lat, lon
success, pnt = tr.TransformPoint(0, col, row, z)
print(pnt)
#  (19.00054957705992, 59.90823713495448, 15.0)
``````