7

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:

enter image description here


Check

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

enter image description here

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:

enter image description here

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?

6
  • The Sentinel image is in an UTM projection - so to plot your point you need to project your point coordinates from LAT/LON to UTM
    – Val
    Sep 2 '21 at 7:49
  • Looking into it, thanks!
    – zabop
    Sep 2 '21 at 7:53
  • (I realized that since I wrote this question, there are more sentinel images available, so if you run the above script now, you might get a different image -> the cross might be at a different place. when using the coordinates I manually set (x=5200 and y=5000). I hope the question is still clear though.)
    – zabop
    Sep 2 '21 at 8:22
  • Once you have the coordinates in map units (so in the UTM projection of the image) you can then convert them to the pixel indices of the array. So your full workflow is lat/lon -> utm -> pixel index
    – Val
    Sep 2 '21 at 8:35
  • @Val, I have an issue with the lat/lon -> utm path, I explaned it here. Looking at JonasV's answer now, thanks Jonas...
    – zabop
    Sep 2 '21 at 8:56
7

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='+')
4
  • Very cool. I think this solved the issue. My crs (for this specific tiff) is 32618, not 32619 - I guess we have just used a different tiff, is this a cause for concern? Visual checking (ie plotting the image and zooming in) seem to suggest that the point I get using your answer is as least as accurate as I can check.
    – zabop
    Sep 2 '21 at 9:16
  • No, it doesn't really matter much. The border for Zone 18N and 19N goes right through -72°E, so in this area both aren't wrong
    – JonasV
    Sep 2 '21 at 9:23
  • Thank you! Would anything change if I want to plot a completely different area? Let's say Tasmania. I'm checking if I can plot a point near Tasmania as well, a few moments and I'm back...
    – zabop
    Sep 2 '21 at 9:24
  • As long as you get the crs and transform from the tiff you are plotting everything should work just fine
    – JonasV
    Sep 2 '21 at 9:26

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.