5

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?

5
  • 2
    What coordinate system is your projected coordinates in? What datum/coordinate system are you lat/lon coordinates being checked against? – Son of a Beach May 3 at 23:10
  • @SonofaBeach I'm really sorry for not being able to explain this much further but this is my very first time dealing with GIS. I know that the photo has such info relating to coordinate systems: it's using WGS 84 and the UTM zone 23S. The EPSG is 32723. I can't really find information for the PCS. Does it help to inform that the closer I'm to the center of the resized image, the more accurate the conversion gets(although it's still off by some relatively big distances.) – PCosme May 3 at 23:48
  • 1
    So the projected coordinate system is WGS 84 / UTM zone 23S as per: epsg.io/32723 . I would guess that your geographic coordinate system against which you're checking your lat/lon coordinates is likely WGS84, but you should confirm that this is the case. – Son of a Beach May 3 at 23:52
  • @SonofaBeach if the GCS that I'm checking against is indeed WGS 84 should the first method I used give me more accurate results? – PCosme May 3 at 23:55
  • 2
    Sorry, I can't comment on that without a whole lot of research that I don't have time for right now. I did once write an iPhone app that could convert between a variety of datums as well as convert between Geographic and UTM projection. But that was a long time ago. These days, if I want to do such conversions, I use GDAL which does it all for me (and it uses PROJ, I believe). These both have Python modules, I think. If you use these, it makes such conversions trivially easy. See: gdal.org and proj.org – Son of a Beach May 4 at 0:07
4

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
3
  • This works partially for what I need, as it does indeed do the proper conversions but it messes the order of my dataframe. My X and Y coordinates are stored in a dataframe, so when I try to store the value of the conversions in it they are all out of place and don't match the respective X and Y values. It returns me ValueError: Length of values (369) does not match length of index (4889) and when I just change from 369 to 4889 it distorts my values. Well. this did lead me to the righ direction though, so thank you a lot. Really. – PCosme May 4 at 16:33
  • 1
    You can try and run the transformation on the entire dataset in one hit if you like: x_hat, y_hat = cell_transform * (df['X'], df['Y']) - this will give you back two Series (x_hat and y_hat) of the same length as your dataframe which you can just slot in as a column. – om_henners May 5 at 23:33
  • Thanks a lot for the help. I actually managed to do it today earlier, that's how I did it (]stackoverflow.com/questions/67403265/…) My solution is in the answer. It's similar to what you described. By the way, would you mind telling me what kind of field should I study to understand this better? Like, how do you know that multiplying the X and Y values by the affine gives me the correct projection? What is the name of this transformation? Affine transformation? – PCosme May 5 at 23:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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