Coordinates doesn't match between metadata and plot

[sorny92]

I just started using EOReader.

I'm following this link: https://eoreader.readthedocs.io/en/v0.19.4/notebooks/windowed_reading.html As I want to visualize a POLYGON over an image of sentinel-2.

The polygon I want to plot is something like this: "POLYGON ((-7.5262454743456715 39.66672692059894, -7.5238403916980525 39.66669637036892, -7.523853372713001 39.666086466637324, -7.52625843422954 39.666117016209796, -7.5262454743456715 39.66672692059894))"

The image I want plot with it has a metadata where the footprint is something like this:

            <Product_Footprint>
                <Global_Footprint>
                    <EXT_POS_LIST>39.744399896499196 -7.832855 39.72445184157785 -6.552124 38.73593643111066 -6.5862427 38.75519833855342 -7.849121 39.744399896499196 -7.832855 </EXT_POS_LIST>
                </Global_Footprint>
            </Product_Footprint>

As you can see my polygon should be contained on the image.

The problem I'm finding is that once I load the image following the code and I plot just the image I get the next coordinates.

As you can see in the previous image the range of x and y coordinates is aroung 1e5-1e6.

My code looks like this:

polygon_str = "POLYGON ((-7.5262454743456715 39.66672692059894, -7.5238403916980525 39.66669637036892, -7.523853372713001 39.666086466637324, -7.52625843422954 39.666117016209796, -7.5262454743456715 39.66672692059894))"
wkt_polygon = wkt.loads(polygon_str)
polygon = gpd.GeoDataFrame(index=[0], crs=prod.crs(), geometry=[wkt_polygon])
polygon = polygon.to_crs(prod.crs())
polygon.plot()

# To load the image I just do this:
prod = Reader().open(filename)
bands = prod.load([GREEN])        
bands[GREEN][:, ::10, ::10].plot()

I'm assuming I'm missing some transformation somewhere as this data is just a .safe downloaded from SNAP.

[remi-braun]

Hey!

I'm glad you want to use EOReader 😄

You maybe don't know, but EOReader is always working with projected (UTM) CRS, to allow computing areas easily. Your Sentinel-2 data is natively in UTM (even if the coordinates you are showing are in lat/lon). So your plot has UTM coordinates 😉

This is a known limitation, and an issue has been adressed here.

---

Moreover, if you are not aware, you cannot use only plot to overlap geographical plots, one will only juxtapose over the other, keeping only the legend of the first graph.

A good way to do this is to use hvplot with geoviews and to do as I've done in this cell.

And don't forget to import this:

import hvplot.pandas
import hvplot.xarray

If you just want to have two plots, just add :

from matplotlib import pyplot as plt

# ... your code before the first plot
plt.show()

# ... your code before the second plot
plt.show()

[sorny92]

Thanks for the fast reply!

I'm just kind of confused because the tile I want to process is on the 30T tile UTM, which lat/lon coordinates match my polygon, but I don't understand what is the range of values of the plot that I attached to the issue.

How do I go around to get the coordinates in UTM once I load my plot? Or maybe what I should do is transform my polygon to the coordinates of my plot, is that something that I should do with EOReader or geopandas should be the way to go?

I assume I'm dropping really simple questions to you but I'm starting in the world of satellite imaging and I'm trying to understand where does my problem come from!

[remi-braun]

No problem for your questions! The coordinates you see are those in UTM which are in meters. If your not familiar with projections, you should read a bit about them and the difference between geographic (lat/lon, i.e. epsg:4326) and projected ones (UTM, i.e. epsg:32[67]\d\d).

If you really want to have lat/lon coordinates in your plot, you should reproject your raster to epsg:4326 before plotting. You can use rioxarray.reproject to achieve it.

[sorny92]

If you really want to have lat/lon coordinates in your plot, you should reproject your raster to epsg:4326 before plotting. You can use rioxarray.reproject to achieve it.

I understand that reproject will actually wrap the whole image so this is not super efficient and I'm better of converting the polygon to UTM, am I right?

Nevertheless, thanks for your help!

[remi-braun]

I understand that reproject will actually wrap the whole image so this is not super efficient and I'm better of converting the polygon to UTM, am I right?

You are already doing it here 😉 :

polygon = polygon.to_crs(prod.crs())

[sorny92]

I understand that reproject will actually wrap the whole image so this is not super efficient and I'm better of converting the polygon to UTM, am I right?

You are already doing it here wink :

polygon = polygon.to_crs(prod.crs())

Well this is part of my confusion. As I understood, this should have work to conver py polygon to the same crs of the image but it keeps it more or less at the same.

I will play with it a bit later.

Thanks for your help!