Geo-coordinates of Sentinel-2 results

Post Reply
martin.boettcher.bc
Posts: 5
Joined: Thu Oct 10, 2019 3:04 pm
company / institution: Brockmann Consult GmbH
Location: Hamburg, Germany

Geo-coordinates of Sentinel-2 results

Post by martin.boettcher.bc »

Dear Francois,

in order to compare results of different processors it is often easiest to merge the outputs of them into one. We do merging with SNAP MergeOp, and we currently work with Sentinel-2 60m resolution. But the Polymer outputs unfortunately refuse to be merged because the pixel positions are slightly different from those SNAP provides for e.g. the C2RCC outputs written to NetCDF. Lat and Lon rasters of Polymer seem to be stretched by the distance of half a pixel into each direction, making the image slightly larger than it should be.

Comparing two products there are always two, and it may be that the other one is wrong. I do not exclude that SNAP is wrong in how it reads the geo-information of the Sentinel-2 L1C product. The reasons that support the SNAP way of reading Sentinel-2 are: SNAP expects that the lat and lon coordinates are centre coordinates. I think that Polymer uses the same assumption because else I would expect that the Polymer image is only shifted in relation to the SNAP output image, but not stretched. Second, I assume that in the L1C product the corner coordinate given in UTM define the north-east corner (!) of the north-east pixel because the coordinates are the same for the resolutions 10m, 20m, and 60m (and I assume the first 6*6 10m pixels cover the same area as the first 60m pixel). This matches with the SNAP way of reading the product and writing 60m processing results.

It seems that Polymer uses the north-east corner coordinates and the calculated south-west corner coordinates of the complete L1C image (by adding 1830*60m) and divides it by 1829 to get the pixel spacing for the "60m" result (that no longer is 60m then). This assumes that the image corner coordinates are pixel centres of the outermost pixels. I do not think that this fits well with the image extend being the same in 10m, 20m and 60m and how I assume the 10m pixel blocks fit into the 60m pixels. Would it be possible to divide by 1830 instead to really get 60m pixels? And would it be possible to shift by half a pixel distance to get pixel centre coordinates in the lat and lon variables. This would make merging with SNAP products possible, and it would make matchups more correct (if my observations are correct).

Best regards,
Martin
User avatar
fsteinmetz
Site Admin
Posts: 306
Joined: Fri Sep 07, 2018 1:34 pm
company / institution: Hygeos
Location: Lille, France
Contact:

Re: Geo-coordinates of Sentinel-2 results

Post by fsteinmetz »

Hi Martin,

Thanks for the detailed feedback !
You are right, it is clear Polymer is wrong in how it reads the lats/lons. As you say, the coordinates of the centre of the north east pixel should not be the same at different resolutions.

Looking at the code (method init_latlon in level1_msi.py), it is clear that an offset should be added. However, my current understanding is that it should be a simple offset without stretching.

Each pixels is currently located in metrics coordinates by :

Code: Select all

    #level1_msi.py:201
    X, Y = np.meshgrid(ULX + XDIM*np.arange(self.totalheight), 
                       ULY + YDIM*np.arange(self.totalwidth))
Where XDIM and YDIM are the sizes of the pixels in metres at the target resolution ; self.totalheight and self.totalwidth are equal to 1830 at 60m resolution. And then these coordinates are projected to geographic coordinates.
To fix the offset, I would add 0.5 to each np.arange.

So, at 60m, the north east pixel would have coordinates (ULX + XDIM*0.5, ULY+YDIM*0.5), and the south west pixel would have coordinates (ULX + XDIM*1829.5, ULY+YDIM*1829.5).

Do you think that would work? There is no division by the grid size, so I can't see where the stretching would come from.

Cheers,
François
martin.boettcher.bc
Posts: 5
Joined: Thu Oct 10, 2019 3:04 pm
company / institution: Brockmann Consult GmbH
Location: Hamburg, Germany

Re: Geo-coordinates of Sentinel-2 results

Post by martin.boettcher.bc »

Dear Francois,

thank you for your quick response. You are right, it is just a shift. I got confused by too many corners I looked at. The fix you propose solves the problem. I have done the following modification to stay with integer calculations, valid since the DIM values are even:

X, Y = np.meshgrid(ULX + XDIM/2 + XDIM*np.arange(self.totalheight),
ULY + YDIM/2 + YDIM*np.arange(self.totalwidth))

Best wishes
Martin
User avatar
fsteinmetz
Site Admin
Posts: 306
Joined: Fri Sep 07, 2018 1:34 pm
company / institution: Hygeos
Location: Lille, France
Contact:

Re: Geo-coordinates of Sentinel-2 results

Post by fsteinmetz »

Great! Good to read that it solves the issue.
I will integrate this modification in further Polymer releases.
Just a minor comment: in python3 you would have to use XDIM//2 to stay with integers, XDIM/2 gives a float. But I don't expect that it should affect the results - maybe the performance ? But very slightly then.
Thanks again!
Cheers,
François
User avatar
fsteinmetz
Site Admin
Posts: 306
Joined: Fri Sep 07, 2018 1:34 pm
company / institution: Hygeos
Location: Lille, France
Contact:

Re: Geo-coordinates of Sentinel-2 results

Post by fsteinmetz »

So, for the record, the patch is:

Code: Select all

diff --git a/polymer/level1_msi.py b/polymer/level1_msi.py
index e2a7b63..c3e1926 100644
--- a/polymer/level1_msi.py
+++ b/polymer/level1_msi.py
@@ -198,8 +198,8 @@ class Level1_MSI(Level1_base):
                 XDIM = int(e.find('XDIM').text)
                 YDIM = int(e.find('YDIM').text)
 
-        X, Y = np.meshgrid(ULX + XDIM*np.arange(self.totalheight), 
-                           ULY + YDIM*np.arange(self.totalwidth))
+        X, Y = np.meshgrid(ULX + XDIM//2 + XDIM*np.arange(self.totalheight), 
+                           ULY + YDIM//2 + YDIM*np.arange(self.totalwidth))
 
         self.lon, self.lat = (proj(X, Y, inverse=True))
yifanaaa
Posts: 3
Joined: Fri Aug 14, 2020 8:59 am
company / institution: Institute of oceanology, Chinese academy of sciences
Location: qingdao, China

Re: Geo-coordinates of Sentinel-2 results

Post by yifanaaa »

So sorry to trouble you. I am using SNAP software to display the corrected msi netcdf data compared with the orginal L1 image.Since I have already modified the level1 msi script as you have mentioned above, shown below. I am really confused of this output.
Image
But the coordinates can not be corresponded. The following image is the relative position of these two images.
Image
You do not have the required permissions to view the files attached to this post.
Post Reply