Sentinel-2 geographic to image coordinates issue

Optical Toolbox
#1

Hi,

I have trouble converting geographic coordinates to image coordinates. My procedure works with most Sentinel-2 products I’ve tried. However, every once in a while, it fails. I’m using Python, but not snappy.

The procedure is as follows.

lat, lon coordinates -> utm coordinates -> inverse projection matrix multiplication -> image coordinates

This is an example that reproduces the issue, final X image coordinate is wrong:

Product data:

Product id: 05303463-cdab-4b95-bc65-5ebaf435d45e
Product name: S2B_MSIL2A_20201112T034009_N0214_R061_T49SBB_20201113T160805
Footprint:
"POLYGON((108.85214171500951 37.35213890138582,108.84875074421133 37.3419731357291,108.79949598951012 37.19538606239605,108.75143165356803 37.04857431885545,108.71390883771207 36.93518822112508,107.63258699686288 36.910119161555244,107.58802752132053 37.8980997878049,108.83567015867037 37.92768043331917,108.85214171500951 37.35213890138582))

Matrix for band B2,B3,B4 (resolution 10). Double checked with product metadata in /GRANULE/L2A_T49SBB_A019249_20201112T034005/MTD_TL.xml and SNAP-generated geotiff using rasterio:

with rio.open('full_image.tif') as dataset:    
    matrix = dataset.transform 

| 10.00, 0.00, 199980.00|
| 0.00,-10.00, 4200000.00|
| 0.00, 0.00, 1.00|

Input data:

- lat, lon coordinates: (37.44, 107.64), somewhere in China.
- utm coordinates: 733558.139047, 4146957.284352

UTM conversion checked with Lat Long to UTM Converter

- inverse matrix:
| 0.10, 0.00,-19998.00|
| 0.00,-0.10, 420000.00|
| 0.00, 0.00, 1.00|
- pixel coordinates:  x= 53357.81, y= 5304.271

Y coordinate seems about right, snap places it about halfway the image height. However, X is five times the image width (10980 pixels).

This is a minimal example code (Python 3.6)

from affine import Affine
import utm

matrix = Affine(10.00,0.0,199980.00,0.0,-10,4200000.00)
print("{}\n".format(matrix))
print("{}\n".format(~matrix))
    
lat = 37.44
lon = 107.64

UTMx, UTMy, zone, band = utm.from_latlon(lat,lon)
print("UTMx = {}, UTMy = {}".format(str(UTMx),str(UTMy)))

x, y = ~matrix*(UTMx, UTMy)
print("x = {}, y = {}".format(str(x),str(y)))

As I said, this code works most of the time, it just seems like the matrix in some products is not right, which can’t be the case, as SNAP/GPT seems to get it right

Thank you

#2

Have you tried more widely used python libraries such as pyproj? You should report your issue upstream .

#3

Well, it seems this is GIS 101, but I just found out. Geospatial imagery is usually referenced to a sinlge utm zone/band, even though it might span across multiple zones. So, a whole product might be referenced to 49N, even though some part of its footprint is on zone 48.

If you perform an automatic coordinate conversion from latlon to utm (i.e: with automatic estimation of UTM zone), a subsequent UTM to pixel coordinates conversion will fail if actual zone of location and the zone product is referenced to do not match.

Long story short:

This conversión from latlon to UTM is correct:

lat = 37.8981
lon = 107.588028

UTMx, UTMy, zone, band = utm.from_latlon(lat,lon)
print("UTMx = {}, UTMy = {}".format(str(UTMx),str(UTMy)))
print("zone: {}, band: {}".format(str(zone), band))

UTMx = 727555.1639196563, UTMy = 4197667.171656173
zone: 48, band: S

But if you want to convert to pixel coordinates afterwards, you must force usage of the utm zone/band product is referenced to:

lat = 37.8981
lon = 107.588028

UTMx, UTMy, zone, band = utm.from_latlon(lat,lon, 49, 'N')
print("UTMx = {}, UTMy = {}".format(str(UTMx),str(UTMy)))
print("zone: {}, band: {}".format(str(zone), band))

UTMx = 199980.0428529873, UTMy = 4200000.022132119
zone: 49, band: N
2 Likes