Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Define MGRS tiles for polar regions #38

Open
johntruckenbrodt opened this issue Jun 22, 2022 · 4 comments
Open

Define MGRS tiles for polar regions #38

johntruckenbrodt opened this issue Jun 22, 2022 · 4 comments
Labels
enhancement New feature or request

Comments

@johntruckenbrodt
Copy link
Member

Currently the processor makes use of a KML file provided for the Sentinel-2 mission to get the characteristics for individual MGRS tiles. This file does not define areas at the poles:

The MGRS itself grid does define tiles in these areas. An alternative to the KML file has to be found to enable processing of scenes in these regions.

@johntruckenbrodt johntruckenbrodt added the enhancement New feature or request label Jul 8, 2022
@johntruckenbrodt
Copy link
Member Author

The MGRS coordinates can be computed using the mgrs library, which would already be an improvement over using the Sentinel-2 KML file:

>>> import mgrs
>>> m = mgrs.MGRS()
>>> print(m.toMGRS(latitude=51, longitude=11.5, MGRSPrecision=0))
'32UPB'
>>> print(m.MGRSToUTM('32UPB'))
(32, 'N', 600000.0, 5600000.0)

UPS polar tile names can also be retrieved using toMGRS but an equivalent to MGRSToUTM for UPS does not exist.

@johntruckenbrodt
Copy link
Member Author

johntruckenbrodt commented Sep 6, 2022

The following can be used to get the geometry of the corresponding tile for a point:

import mgrs
from spatialist.vector import bbox

lat = 84
lon = 11.5

m = mgrs.MGRS()
mgrs_tile = m.toMGRS(latitude=lat, longitude=lon, MGRSPrecision=0)
zone, hem, east, north = m.MGRSToUTM(mgrs_tile)

epsg = 32600 + zone
if hem == 'S':
    epsg += 100

coordinates = {'xmin': east, 'xmax': east + 109800,
               'ymin': north + 100000, 'ymax': north - 9800}

with bbox(coordinates=coordinates, crs=epsg) as geom:
    print(geom.convert2wkt(set3D=False)[0])

This prints the geometry for tile 33XVP:

POLYGON ((400000 9400000,400000 9290200,509800 9290200,509800 9400000,400000 9400000))

which is slightly different to the geometry in the Sentinel-2 KML file:

POLYGON ((399960 9400020,399960 9290220,509760 9290220,509760 9400020,399960 9400020))

However, not all tiles show this difference.

Also, in areas of UTM zone overlap not all MGRS tiles are covered by the S2 KML. This is further described here:
https://gis.stackexchange.com/questions/227847/finding-sentinel-tile-for-specific-long-lat-coordinate

@johntruckenbrodt
Copy link
Member Author

See here for an explanation of the point above: opendatacube/odc-stac#94 (comment)

@johntruckenbrodt
Copy link
Member Author

Here is a first draft:
image

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request
Projects
None yet
Development

No branches or pull requests

1 participant