How to find a particular coordinate on a Sentinel image using proj:transform?

Question

Background

I am retrieving images from the Sentinel satellites. I plot these images, using pixel coordinates. I would like to show points specified using latitude and longitude on those images.

---

Data acquisition

I am using pystac-client to retrieve Sentinel images. After installing the module (pip install pystac-client), I do, following the readme:

from pystac_client import Client
catalog = Client.open("https://earth-search.aws.element84.com/v0")

mysearch = catalog.search(collections=['sentinel-s2-l2a-cogs'],
                          bbox=[-72.5,40.5,-72,41],
                          query =  {"eo:cloud_cover":{"lt":1}},
                          max_items=10)   

print(f"{mysearch.matched()} items found")

I save the result as a dict by:

resdict = mysearch.get_all_items_as_dict()

For example, select a single result:

resdict['features'][3]['assets']['B08']

This above is a dict:

{'eo:bands': [{'center_wavelength': 0.8351,
   'common_name': 'nir',
   'full_width_half_max': 0.145,
   'name': 'B08'}],
 'gsd': 10,
 'href': 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/19/T/BF/2021/8/S2A_19TBF_20210827_0_L2A/B08.tif',
 'proj:shape': [10980, 10980],
 'proj:transform': [10, 0, 199980, 0, -10, 4600020, 0, 0, 1],
 'roles': ['data'],
 'title': 'Band 8 (nir)',
 'type': 'image/tiff; application=geotiff; profile=cloud-optimized'}

---

Plot result

Imports for plotting:

import rasterio
import matplotlib.pyplot as plt

Plot the .tif file mentioned in the above dict (the href entry there):

plt.figure(figsize=(8,8))
src = rasterio.open(resdict['features'][3]['assets']['B08']['href'])
plt.imshow(src.read(1), cmap='pink')

Result:

---

Check

Check if this is indeed what's expected with Google Maps (using coordinates from the bbox kwarg in the catalog.search() function above):

The same geographical area appears, so what we are plotting is indeed the tip of Long Island, just as we wanted (based on bbox values).

---

Problem

I would like to be able to do what Google Maps did above: take a point's latitude and longitude, and plot it on my image (ie obtain this point's pixel coordinates on my satellite image).

---

What I know:

The dictionary above, which contained the url for the .tif file, has the proj:transform entry. In our case, it is:

[60, 0, 199980, 0, -60, 4500000, 0, 0, 1]

These 9 numbers are mentioned in this issue:

... 9 elements in the array (3x3 matrix) because it's accounting for 3 dimensions, but the last row is just "0 0 1". rasterio includes it when getting the transform via "src.transform"

So the numbers which matter seem to be the first 6. This is reassuring, since affine, a package "describing affine transformation of the plane" requires 6 numbers for such transformation: a, b, c, d, e, f. I suspect that the first 6 numbers in proj:transform are these 6 numbers.

_{Tangential to the topic & I'm just learning yet about projections, but it is probably no accident that there are 2^2*(2+1)/2 = 6 independent Christoffel coefficients for a 2-dimensional metric.}

---

What I don't know

I don't know how to use those 6 numbers (ie 60, 0, 199980, 0, -60, 4500000) to plot like Google Maps does. Given the input 41,-72.5 (the coordinates I gave to Google Maps, see above screenshot), I would like to get something like this:

Where the cross is exactly in the same geographical location as the red pointer on Google's map.

(This figure above was produced with this script:

plt.figure(figsize=(8,8))
src = rasterio.open(resdict['features'][3]['assets']['B08']['href'])
plt.imshow(src.read(1), cmap='pink')

# BELOW IS JUST AN APPROXIMATION
x = 5200 # HOW TO GET THESE PIXEL COORDINATES EXACTLY?
y = 5000 # HOW TO GET THESE PIXEL COORDINATES EXACTLY?

plt.scatter([x],[y],s=10000,c='r',marker='+')

x and y are just visual approximations.)

---

The question

Given a .tif image and the numbers in proj:transform as described above, how can I map latitude-longitude coordinates to pixel coordinates on my image?

Answer 1

So what you basically want to do is transform the coordinates you give in one coordinate reference system (CRS) to another CRS.

With rasterio the CRS of a GeoTiff is saved in the attribute crs which you can easily access:

src = rasterio.open(resdict['href'])
src.crs

which returns CRS.from_epsg(32619). For Google Maps you need to know that the CRS is WGS84, which has the code EPSG:4326 (your ordinary lat/lon).

Now we have everything we need to transform from one CRS to another:

from rasterio.crs import CRS
from rasterio import warp

dst_crs = src.crs
src_crs = CRS.from_epsg(4326)

lat = 41.15
lon = -71.6

x, y = warp.transform(src_crs, dst_crs, [lon], [lat])
print(x, y)

The last step is transforming this to pixel indices which we can use to plot:

from rasterio import transform

row, col = transform.rowcol(src.transform, x, y)
print(row, col)

plt.imshow(src.read(1), cmap='pink')

plt.scatter([col],[row],s=10000,c='r',marker='+')