Wrong tile for a specific point

[tonyVeco]

I have used the mgrs library to get the tile for this point :

Lat = 39.9770 Lon = 116.3810

I get this tile:

50SMK4714425387

this does not correspond to any exsisting Sentinel2 tile. It should be 50TMK

[DanielJDufour]

Hi, @tonyVeco . Thank you for submitting this issue! I definitely want to fix this issue ASAP!

This issue happens when a 100km grid square crosses a latitude boundary and some of its squares are in latitude band S (32° N to 40° N) and some are in latitude band T (40° N to 48° N). The getLetterDesignator function used by this library returns the latitude band in which the point resides. For a better explanation of this issue read: https://en.wikipedia.org/wiki/Military_Grid_Reference_System#Squares_that_cross_a_latitude_band_boundary

However, you rightly point out that this presents challenges for finding Sentinel 2 imagery tiles. I will create a solution for you. First, I will have to better understand how sentinel assigns a latitude band to a 100km grid zone. Is it based on the South-West point or center of the 100km grid square? Additionally, if I understand you correctly, you are just using this library to get the 100km grid square id and don't care about the rest after it (e.g. 4714425387)?

Thank you again for bringing this to our attention. I'd love to make this work for you!

[tonyVeco]

Dear Daniel,

unfortunately I don't know how it works for the latitude band of the Sentinel2 tiles, I am really sorry I can not help you with that. Yes actually I need only the tile ID because I have to search into a product catalogue that is based on the Tile ID, I don't need to go into internal coordinates . Having this problem solved would save me because I am using it for a project and I really need it. Thanks for helping!!!

[tonyVeco]

Dear Daniel,

I am really appreciating your help inn solving this unexpected issue. Do you think it could be solved in the next few weeks? your mgrs library will be used in some Jupyter Notebooks for a current international project. Also let me know how we could credit you in these notebooks.

Regards

Antonio

[DanielJDufour]

@tonyVeco , this will definitely be fixed in the next few weeks. No credit required :-)

[Klaus-Tockloth]

The very matured GeoConvert (1.49) also results in "50SMK4714425387" ...

$ echo 39.9770, 116.3810 | GeoConvert -m
50SMK4714425387

Where came the information from, that the result should be "50TMK"?

[tonyVeco]

Actually I had an error message from a web service where I submitted a query specifying both the coordinates of the point and the tile ID 50SMK, then I tried another service to look directly the Sentinel2 tiles and actually the 50SMK seems to be not the correct one

[Klaus-Tockloth]

The red line in the picture above shows that the 10km square is divided into two part. Coords in the part above the red line are starting with 50TMK (Cross = Center 50TMK47132871). Coords in the part below the red line are starting with 50SMK (Marker = Cursor = 50SMK47162556). The given coords "39.9770 116.3810" (the blue marker) are within the part below. What you expect for the Sentinel2 IDs is, that each coord in the 10km square has to start with "50TMK". Right? Is there an algorithm where Sentinel2 describes how one has to handle this special case?

[tonyVeco]

this is exactly the point said by Daniel, I think he will fix it :) I don't know the algorithm about Sentinel2 tiling..

[Klaus-Tockloth]

Just for clarification: I currently don't see any problem, the library is working correctly.

[tonyVeco]

The library is working correctly but it has not a 100% correspondance with the Sentinel2 tiling grid, so I asked it if it possible to do something for it. For clarification, download this : https://sentinel.esa.int/documents/247904/1955685/S2A_OPER_GIP_TILPAR_MPC__20151209T095117_V20150622T000000_21000101T000000_B00.kml it is the S2 tiling grid and you will not find the 50SMK

[Klaus-Tockloth]

That's what the referenced kml file shows:

Sentinel2 provides probably only full square images. That makes sense. But the MGRS derived image/tile name "50TMK" is only correct for a part of the image ... and that's your problem. The center of Beijing has a MGRS that starts with "50SMK".

[tonyVeco]

exactly, the Tiling grid of S2 does not contain any 50SMK Id, and this is the problem. In some way some specific modification seems to be needed with the MGRS package to work with the few exceptions in the S2 tiles… Any other solution?

[Klaus-Tockloth]

Any other solution?

Build and use your own lookup table. Everything you need is in the referenced kml file. For given coords the lookup should then result in:

• no image found (poles)
• 1 image found (coords are within the core of an image)
• 2 images found (coords are on an overlapping margin)
• 4 images found (coords are on an overlapping corner)

Data example for Beijing (WKT format): TILE_ID = 50TMK MULTIPOLYGON(((115.816827393401 40.644794803776,115.83386830445 39.655752572233,117.113779914526 39.6615512269476,117.115443067619 40.6507988109166,115.816827393401 40.644794803776)))

[DanielJDufour]

Hi, @tonyVeco and @Klaus-Tockloth . Thank you for all your research. I did a little more research and wanted to add this to the discussion. I downloaded the 100km Square Identifier Polygons from NGA at http://earth-info.nga.mil/GandG/update/coordsys/data/mgrs/world_mgrs.zip and loaded it into QGIS. NGA cuts the 100km squares into two where they cross a latitude boundary.

Here's a zoomed in picture with the point in question:

My next step will be to write a script to better understand how Sentinel 2 does its tiling. I'd like to understand better how the Sentinel 2 grid is created before offering an opinion. I appreciate your patience with me :-)

[DanielJDufour]

For anyone else interested, here's the relevant section from the KML file mentioned above:

<Placemark>
    <name>50TMK</name>
    <Snippet maxLines="0"></Snippet>
    <description><![CDATA[TILE PROPERTIES<br><table border=0 cellpadding=0 cellspacing=0  width=250 style="FONT-SIZE: 11px; FONT-FAMILY: Verdana, Arial, Helvetica, sans-serif;"><tr><td bgcolor="#E3E1CA" align="right"><font COLOR="#000000"><b>TILE_ID</b></font></td><td bgcolor="#E4E6CA"> <font COLOR="#008000">50TMK</font></td></tr><tr><td bgcolor="#E3E1CA" align="right"><font COLOR="#000000"><b>EPSG</b></font></td><td bgcolor="#E4E6CA"> <font COLOR="#008000">32650</font></td></tr><tr><td bgcolor="#E3E1CA" align="right"><font COLOR="#000000"><b>MGRS_REF</b></font></td><td bgcolor="#E4E6CA"> <font COLOR="#008000">40.64479965 115.81730038 40.650856516 117 39.749907519 117 39.744039563 115.83284812</font></td></tr><tr><td bgcolor="#E3E1CA" align="right"><font COLOR="#000000"><b>UTM_WKT</b></font></td><td bgcolor="#E4E6CA"> <font COLOR="#008000">MULTIPOLYGON(((399960 4500000,399960 4390200,509760 4390200,509760 4500000,399960 4500000)))</font></td></tr><tr><td bgcolor="#E3E1CA" align="right"><font COLOR="#000000"><b>LL_WKT</b></font></td><td bgcolor="#E4E6CA"> <font COLOR="#008000">MULTIPOLYGON(((115.816827393401 40.644794803776,115.83386830445 39.655752572233,117.113779914526 39.6615512269476,117.115443067619 40.6507988109166,115.816827393401 40.644794803776)))</font></td></tr></table>]]></description>
    <styleUrl>#MapRollOver</styleUrl>
    <MultiGeometry>
        <Polygon>
            <outerBoundaryIs>
                <LinearRing>
                    <coordinates>
                        115.8168273934,40.6447948038,0 117.1154430676,40.6507988109,0 117.1137799145,39.6615512269,0 115.8338683045,39.6557525722,0 115.8168273934,40.6447948038,0 
                    </coordinates>
                </LinearRing>
            </outerBoundaryIs>
        </Polygon>
        <Point>
            <coordinates>116.4661352305,40.1532756916,0</coordinates>
        </Point>
    </MultiGeometry>
</Placemark>

[DanielJDufour]

If it's helpful for new people catching up with this issue, here's another visualization, with all the relevant data on it:

[DanielJDufour]

Here's an example of a tile in the KML file that is assigned to the area underneath it's latitude band boundary where most of the tile resides:

<Placemark>
    <name>50NMP</name>
    <Snippet maxLines="0"></Snippet>
    <description><![CDATA[TILE PROPERTIES<br><table border=0 cellpadding=0 cellspacing=0  width=250 style="FONT-SIZE: 11px; FONT-FAMILY: Verdana, Arial, Helvetica, sans-serif;"><tr><td bgcolor="#E3E1CA" align="right"><font COLOR="#000000"><b>TILE_ID</b></font></td><td bgcolor="#E4E6CA"> <font COLOR="#008000">50NMP</font></td></tr><tr><td bgcolor="#E3E1CA" align="right"><font COLOR="#000000"><b>EPSG</b></font></td><td bgcolor="#E4E6CA"> <font COLOR="#008000">32650</font></td></tr><tr><td bgcolor="#E3E1CA" align="right"><font COLOR="#000000"><b>MGRS_REF</b></font></td><td bgcolor="#E4E6CA"> <font COLOR="#008000">8.1410194048 116.09227432 8.142034146 117 7.2374656175 117 7.2365648459 116.09419428</font></td></tr><tr><td bgcolor="#E3E1CA" align="right"><font COLOR="#000000"><b>UTM_WKT</b></font></td><td bgcolor="#E4E6CA"> <font COLOR="#008000">MULTIPOLYGON(((399960 900000,399960 790200,509760 790200,509760 900000,399960 900000)))</font></td></tr><tr><td bgcolor="#E3E1CA" align="right"><font COLOR="#000000"><b>LL_WKT</b></font></td><td bgcolor="#E4E6CA"> <font COLOR="#008000">MULTIPOLYGON(((116.09191125891 8.14101859288788,116.09400773071 7.14792556358065,117.08839321863 7.14880742694774,117.088597794983 8.14202447886391,116.09191125891 8.14101859288788)))</font></td></tr></table>]]></description>
    <styleUrl>#MapRollOver</styleUrl>
    <MultiGeometry>
        <Polygon>
            <outerBoundaryIs>
                <LinearRing>
                    <coordinates>
                        116.0919112589,8.1410185929,0 117.088597795,8.1420244789,0 117.0883932186,7.1488074269,0 116.0940077307,7.1479255636,0 116.0919112589,8.1410185929,0 
                    </coordinates>
                </LinearRing>
            </outerBoundaryIs>
        </Polygon>
        <Point>
            <coordinates>116.5902545269,7.6449750212,0</coordinates>
        </Point>
    </MultiGeometry>
</Placemark>

[DanielJDufour]

Here's a picture taken from QGIS:

[DanielJDufour]

This information leads me to believe that ESA assigns the latitude band to a tile based on the centroid of the tile. This won't be confirmed until we write a script to test this hypothesis.

I'm going to experiment and see how hard it would be to incorporate an ability to specify a precision level of 100km when running forward and resolving latitude band based on the centroid. forward currently only accepts an accuracy argument to 10000 meters.

Please let me know if I have misunderstood anything.

[DanielJDufour]

@tonyVeco and @Klaus-Tockloth , please consult this pull request and let me know if it makes sense to you.

I'd like to add support for 100km grid square resolution using the centroid without impacting the code or usability much and I think I have found a common ground.

Looking forward to your thoughts :-)

[tonyVeco]

Hi Daniel,

actually I am not an expert on these kind of topics , I work with geospatial data, but more on the processing side than on the geometric complexity behind them. What I could in the next 3 days, is to ask for info to some Sentinel2 expert, I am attending the ESA Living Planet Symposium and I think I can get some expert to have clarification about it. It seems to me reasonable the centroid option, and I would be happy to test the modified mgrs package.. But I will also try to get more info for you about it. Thanks in advance Antonio

[DanielJDufour]

Per the discussion on this pull request, I'm closing this issue. Unfortunately, the S2 tilling grid's latitude band lettering does not conform to a consistent approach, so this library would not be a useful solution to this problem.

Thank you for the discussion and helping identify some interesting things about Sentinel 2 tiling. I'm hoping we see another solution provided outside this library! :-)

[muety]

Super interesting discussion here! Thanks a lot for all your research, @DanielJDufour, @Klaus-Tockloth. I just came across the exact same issue, so this was very helpful to me.

Conclusion is that the only way to reliably find the S2 image for a given lat/lon is to build a lookup table from the above KML file (and this way more or less "ditch" the MGRS concept)? Is that correct? If yes, has anyone already come up with a library for that?

Also, the above link to the MGRS 100km² world grid seems to be broken. Does anyone still have that grid stored somewhere? I'd really like to have it for further experimenting and better understanding the issue.

[DanielJDufour]

Hi, @muety . Yes, that's correct. I think this is the correct link now for what you are looking for: https://earth-info.nga.mil/php/download.php?file=coord-mgrs. You can find it on this page under the Coordinate Systems -> Data tab: https://earth-info.nga.mil/index.php

Unfortunately, I'm not aware of any JavaScript library focused on getting the S2 info correct for a given lat/lon. The way I've seen most people solve this is loading up the KML in Python and doing a point in polygon calculations there. But please know, I'm not aware of all the S2 work out there, so perhaps someone else has solved this.

A sentinel 2 library for this would be a great open source contribution if you wanted to work on it :-)

Let me know if you have any other questions.

[DanielJDufour]

Lastly, I would check out how Sentinel Hub is doing things on their Github. They are usually leading the way on this sort of thing.