Translating X, Y coordinates to lat/long

Question

I have X and Y coordinates from a resized image that has its origins from a Sentinel 2B photo. The original photo is 109.681 by 110.334 km. The re-dimensioned file is 370 by 369. I also know the latitude and longitude coordinates from the upper right and left corner, center and lower right and left corner.

So by multiplying a desired X and Y pixel by 300m(how much 1 pixel equals to) and applying the Euclidean distance formula I can find out how far said pixel is from my reference point(I'm using the upper right corner as reference.) My first attempt at translating the Cartesian coordinates to lat and long went as following:

def pixel_meter():
    global x_meter
    global y_meter
    x_pixel = int(input('Which X pixel? '))
    y_pixel = int(input('Which Y pixel? '))
    x_meter = (x_pixel * PIXEL_METER) * -1
    y_meter = (y_pixel * PIXEL_METER) * -1


def conversion():
    dLat = y_meter / EARTH_RADIUS

    dLon = x_meter / (EARTH_RADIUS * np.cos(np.pi * LAT_REFERENTIAL / 180))

    lat_final = LAT_REFERENTIAL + dLat * 180 / np.pi
    lon_final = LON_REFERENTIAL + dLon * 180 / np.pi

    print(lat_final, lon_final)

That got me close, but when I checked it on a map online it was still off by some kilometers to what would make sense by judging from the pixels displayed by the resized picture.

Then after searching some more I found about the Vincenty's formula and how more precise it is. With pygeodesy I tried this:

from pygeodesy.ellipsoidalVincenty import LatLon


p = LatLon(-18.088506124573, -44.907758947306) #Reference point
d = p.destination(85854, 224.47) # Distance to my desired X, Y coordinates and bearing from center.
print (d.lon, d.lat)

The distance is calculated using Euclidean's distance.

The result is still off, and I'm certain it is because of the bearing. I don't know, and found no way to know what would be the bearing of my points of interest(more than 4 thousand X,Y points).

Is there any way to calculate the bearing just from the information I have?

That is, latitude and longitude from center, upper corners, lower corners and original covered area. I'm open to any possible methods, really.

Also, if it's not possible to calculate the bearings then is there any other way to make this translation between coordinates more precise?

Answer 1

Based on your comments there are two steps here:

1. Translating your known coordinates in latitude and longitude to the projected coordinates of EPSG:32723.
2. Calculating the distance to a given cell.

To start with you'll need the folowing Python libraries installed: pyproj, affine and numpy. There are a lot of other Python libraries for working with spatial data you should definitely check out (GDAL, rasterio, fiona, geopandas, etc.), but I'll go for a quick solution.

First, start by defining some projections:

import pyproj

utm23s = pyproj.CRS('epsg:32723')
wgs84 = pyproj.CRS('epsg:4326')  # the EPSG code for WGS84 (lat/lng) is 4326

transformer = pyproj.Transformer.from_crs(wgs84, utm23s, always_xy=True)

I'll leave you to look up how projections work, but the easiest way to think about it is a math to transform positions on the globe to a flat plane.

Next step is to calculate the coordinate of the upper left coordinate array.

lat, lng = -18.088506124573, -44.907758947306
ul_x, ul_y = transformer.transform(lng, lat)

This can be used to create an affine object (a "2d linear mapping") that is used to transform cell coordinates to projected coordinates in metres.

import affine

cell_transform = affine.Affine(
    300,  # w-e pixel resolution / pixel width.
    0,  # row rotation (typically zero).
    ul_x,  # x-coordinate of the upper-left corner of the upper-left pixel.
    0,  # column rotation (typically zero).
    -300,  # n-s pixel resolution / pixel height (negative value for a north-up image).
    ul_y  # y-coordinate of the upper-left corner of the upper-left pixel.
)

From here you can transform from cell coordinates to projected coordinates and vice versa easily.

projected_x, projected_y = cell_transform * (column, row)
column, row = ~cell_transform * (projected_x, projected_y)

You can transform your cells of interest and use simple Euclidian distance to each location one at a time, or each cell in bulk if you prefer:

xx, yy = cell_transform * np.meshgrid(np.arange(370), np.arange(369))

distance = ((xx - ul_x)**2 + (yy - ul_y)**2)**0.5