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

Playing around with dynamic specifications from wkt2 strings #2

Merged
merged 6 commits into from
Dec 13, 2022

Conversation

jlaura
Copy link
Contributor

@jlaura jlaura commented Dec 7, 2022

Went with urllib to keep everything in the standard library for pulling the input WKT2.

Copy link
Owner

@AndrewAnnex AndrewAnnex left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

oops thought I was committing a change

@AndrewAnnex
Copy link
Owner

AndrewAnnex commented Dec 8, 2022

weird, so my fix apparently broke the tests because sometimes a None is returned for the authority. I am digging into it now but it looks like converting to a dict first will always resolve the id

for example the CRS for Steins (2015) - Sphere / Ocentric has a valid authority of ('IAU_2015', '200286700') but for Toutatis (2015) - Sphere / Ocentric, I get a None for the authority, but when converted to a dict first I get {'authority': 'IAU', 'code': 200417900, 'version': 2015} which seems fine

@AndrewAnnex
Copy link
Owner

weird, so here are the two wkts:

GEOGCRS["Toutatis (2015) - Sphere / Ocentric",\n    DATUM["Toutatis (2015) - Sphere",\n    \tELLIPSOID["Toutatis (2015) - Sphere", 1331.6666666666667, 0,\n\t\tLENGTHUNIT["metre", 1, ID["EPSG", 9001]]]],\n    \tPRIMEM["Reference Meridian", 0,\n            ANGLEUNIT["degree", 0.0174532925199433, ID["EPSG", 9122]]],\n\tCS[ellipsoidal, 2],\n\t    AXIS["geodetic latitude (Lat)", north,\n\t        ORDER[1],\n\t        ANGLEUNIT["degree", 0.0174532925199433]],\n\t    AXIS["geodetic longitude (Lon)", east,\n\t        ORDER[2],\n\t        ANGLEUNIT["degree", 0.0174532925199433]],\n\tID["IAU", 200417900, 2015],\n\tREMARK["Use R_m = (a+b+c)/3 as mean radius. Use mean radius as sphere radius for interoperability. Source of IAU Coordinate systems: doi://10.1007/s10569-017-9805-5"]]

and

GEOGCRS["Steins (2015) - Sphere / Ocentric",\n    DATUM["Steins (2015) - Sphere",\n    \tELLIPSOID["Steins (2015) - Sphere", 2700, 0,\n\t\tLENGTHUNIT["metre", 1, ID["EPSG", 9001]]],\n\t\tANCHOR["Topaz : 0"]],\n    \tPRIMEM["Reference Meridian", 0,\n            ANGLEUNIT["degree", 0.0174532925199433, ID["EPSG", 9122]]],\n\tCS[ellipsoidal, 2],\n\t    AXIS["geodetic latitude (Lat)", north,\n\t        ORDER[1],\n\t        ANGLEUNIT["degree", 0.0174532925199433]],\n\t    AXIS["geodetic longitude (Lon)", east,\n\t        ORDER[2],\n\t        ANGLEUNIT["degree", 0.0174532925199433]],\n\tID["IAU", 200286700, 2015],\n\tREMARK["Use mean radius as sphere radius for interoperability. Source of IAU Coordinate systems: doi://10.1007/s10569-017-9805-5"]]

I see Id's in both of these but Steins has an anchor defined for the ellipsoid while the first (failing) one doesn't

@jlaura
Copy link
Contributor Author

jlaura commented Dec 8, 2022

Those are some legit gymnastics in order to get the codes! Glad you got that figured out.

I want to look at TiTiler / Morecantile and determine how best to associate the input data proj with the identifier. For example IAU_30100_2015 is not going to be meaningful to 99.9% of the users.

I sent the link to the test TiTiler deploy (AWS) via email. At the /docs URL, you can see the swagger spec/UI. In there, my first approach was to add body specific endpoints (which is not ideal for a myriad of reasons, but works for now). So, <link>/europa/tiles/<IAU_#####_2015>/{z}/{x}/{y} would be one endpoint that was usable to tile one's data using the service if we just straight ingested the output of generate.py. Thoughts?

@jlaura
Copy link
Contributor Author

jlaura commented Dec 10, 2022

Messing with this a bit more this morning. Here are a few things that I noticed:

    import morecantile
    tms = morecantile.TileMatrixSet.parse_file('Moon_2000_tms.json')
  • Copying one of these custom TMS *.json files into the morecantile data area, which seems to be one method to register them causes the same error as above.

Option 1: Hack the CRSes in the JSON files.

In the short term, it might be most elegant to perform this hack (replace None with 4326 since these are all GCS. It would be significantly better to figure out who maintains opengis.net and try to get the IAU authorities added.

Edit: The above is partially correct. Randomly browsing opengis.net I see an endpoint for IAU! So, the hack is unnecessary in four cases. What needs to be done is to get the correct supportedCRS and boundingBox:crs links for the code. opengis.net has entries for Mars2000, Mercury2000, Moon2000, and Venus2000 (which is kind of random...).

If one hacks the supportedCRS and and boundingBox:crs it is the validator works. The resulting TMS looks like:

{
    "type": "TileMatrixSetType",
    "title": "Sun (2015) - Sphere / Ocentric",
    "identifier": "IAU_1000_2015",
    "supportedCRS": "http://www.opengis.net/def/crs/EPSG/0/4326",
    "boundingBox": {
        "type": "BoundingBoxType",
        "crs": "http://www.opengis.net/def/crs/EPSG/0/4326",
        "lowerCorner": [
            -90.0,
            -180.0
        ],
        "upperCorner": [
            90.0,
            180.0
        ]
    },

Option 2: Register a custom TMS to morecantile

  • It is also possible to use the following to properly create a custom TMS (note that this is showing examples with the ESRI authority and the IAU authorities):
from pyproj import CRS

mars_crs = CRS.from_authority("ESRI", 104971)
ESRI104971 = morecantile.TileMatrixSet.custom(
    crs=mars_crs, 
    identifier="Mars_2000_Sphere_ESRI",
    extent=mars_crs.area_of_use.bounds, 
    matrix_scale=[1, 1],
)

iau_mars_crs = CRS.from_authority("IAU", "49900")
IAU49900 = morecantile.TileMatrixSet.custom(
    crs=iau_mars_crs, 
    identifier="Mars_2000_Sphere_IAU",
    extent=mars_crs.area_of_use.bounds, 
    matrix_scale=[1, 1],
)
  • I played around with parsing the WKT into custom TMSes for morecantile. This seems to be working to get these registered.
import urllib
with urllib.request.urlopen('https://raw.githubusercontent.com/pdssp/planet_crs_registry/main/data/result.wkts') as response:
    resp = response.read().decode(response.headers.get_content_charset())
    
import os
geocrss = []
for wkt_str in resp.split(os.linesep + os.linesep):
    if 'GEOGCRS' in wkt_str[:7] and 'TRIAXIAL' not in wkt_str:  # insanely hacky
        geocrss.append(wkt_str)
        
for crs in geocrss:
    crs_obj = CRS(crs)
    title = crs_obj.name
    auth = crs_obj.to_authority(min_confidence=25)
    if auth is not None:
        authority_version, code = auth
        authority, version = authority_version.split('_')
        identifier = f'{authority}_{code}_{version}'

        if crs_obj.area_of_use:
            bounds =  crs_obj.area_of_use.bounds  # Not sure that this is working?
        else:
            bounds = [-180, -90, 180, 90]
        custom_tms = morecantile.TileMatrixSet.custom(
                                                    crs=crs_obj, 
                                                    identifier=identifier,
                                                    extent=bounds,
                                                    matrix_scale=[1, 1])
        
        # Strange append because morecantile.tms is a custom object
        morecantile.tms = morecantile.tms.register(custom_tms)
morecantile.tms.list()

Results in (truncated):

['LINZAntarticaMapTilegrid',
 'EuropeanETRS89_LAEAQuad',
 'CanadianNAD83_LCC',
 'UPSArcticWGS84Quad',
 'NZTM2000',
 'NZTM2000Quad',
 'UTM31WGS84Quad',
 'UPSAntarcticWGS84Quad',
 'WorldMercatorWGS84Quad',
 'WGS1984Quad',
 'WorldCRS84Quad',
 'WebMercatorQuad',
 'IAU_1000_2015',
 'IAU_19900_2015',
 'IAU_19901_2015',
 'IAU_29900_2015',
 'IAU_30100_2015',
 'IAU_39900_2015',
 'IAU_39901_2015',
 'IAU_40100_2015',
 'IAU_40200_2015',
 'IAU_49900_2015',
...

All of this is to say: How are we anticipating one use these custom TMS?

@AndrewAnnex
Copy link
Owner

@jlaura That is a lot of great detective work and I can't respond to it all now, but a few things I can say:

  1. Registering the custom TMSs is already a part of planetcantile

    import importlib.resources as importlib_resources
    import pathlib
    from morecantile import TileMatrixSet
    from morecantile import tms
    def to_tms(path):
    with importlib_resources.as_file(path) as src:
    tms = TileMatrixSet.parse_file(src)
    return tms
    planet_tms_jsons = importlib_resources.files('planetcantile.data').glob('*.json')
    planet_tmss = map(to_tms, planet_tms_jsons)
    planetary_tms = tms.register(*planet_tms_jsons, overwrite=True)

    Which was the initial intent of this repo to both generate correct TMS specs as json files and to register them in python for those using python. I hope there are also other projects (JS frontends?) to contribute the json files elsewhere so that other non-python frameworks can use them or we could also make additional libraries to distribute these files. Section Annex D of the TMS 2.0 spec (https://docs.opengeospatial.org/is/17-083r4/17-083r4.html#toc48) lists some common ones for Earth so I can imagine maybe a similar list being published for planetary bodies along with the json files.. My hope is that these custom TMS become standardized and just as easy to access/use as 4326.

  2. That's really interesting about the opengis.net, since it is the OGC definitions server we should do some work to identify what changes they need to support IAU better and see how quickly that can happen. I think I identified another specification (well known scale sets) that need to be adjusted to each planetary body and I think those are also defined at opengis.net (see https://defs.opengis.net/vocprez/object?uri=http%3A//www.opengis.net/def/wkss). The alternative would be to do the 4326 hack, but let's see how far we can go without it.

  3. I planned to add some steps to the github action to add any changes to the json files made in a PR as a commit to the PR also so the generate.py file is only really used by the ci service to support continuous delivery of these specs.

@AndrewAnnex
Copy link
Owner

@jlaura this looks relevant to our discussion opengeospatial/NamingAuthority#212

@jlaura
Copy link
Contributor Author

jlaura commented Dec 11, 2022

Good find! This is what Trent has been talking about. I did not realize the timeline and that it would be linking through opengis.net. That solves the problem on our side as far as I can tell!

We'll have to check morecantile/proj to ensure that they are encoding the URL properly, e.g., that they didn't hard code the /0/ where we need the IAU year /2015/ and IAU vs EPSG. That should get the identifiers straightened out and the code up and running. 🚀

Tangential: Here is where I've been iterating with a STAC proj extension maintainer about getting the codes supported in there. All these small friction points will need either good docs on our side or changes on projects in order to make use of things like the official IAU codes.

@AndrewAnnex
Copy link
Owner

AndrewAnnex commented Dec 11, 2022

as an update to point 2 I made above, I made this PR for morecantile to support generating more-correct scaleDenominators by using the CRS ellipsoid's semi major axis rather than hard-coding the earth radius (developmentseed/morecantile#92). Looking into it the WKSS part may be partially deprecated in TMS 2.0 so we may not need to define additional wkss's

I totally agree about figuring out a way to remove NAIF Ids from the names of the bodies, but I am unsure what the right solution would be for that at the moment if the NAIF id is used is the code and if the code must be an integer?

for the urls encodings in morecantile this looks like the relevant bit of code
https://github.com/developmentseed/morecantile/blob/689cad98914033e746f7eca78719c6e751f843a4/morecantile/models.py#L72-L75

@AndrewAnnex
Copy link
Owner

AndrewAnnex commented Dec 11, 2022

I'll make a pr for that also as it could use the same logic I added to this pr to grab the code/version/authority (edit: done! developmentseed/morecantile#93)

@jlaura
Copy link
Contributor Author

jlaura commented Dec 11, 2022

😆 I was just logging on to get that pushed in. 👍 Glad to see the PRs opened.

With a hopeful update to opengis, I think that these will be all set to deploy without a bunch of hackiness. Once that is the case, I am happy to host these as a lambda at a USGS URL as a community service.

@vincentsarago
Copy link

👋 this is a great conversation 😄
Let me know if I can help with anything. I've seen mentioned to the TMS spec 2.0, just to be clear right now morecantile creates TMS of version 1.0, if it can be helpful to you I can start working of switching to 2.0 🤷‍♂️

@AndrewAnnex
Copy link
Owner

@vincentsarago thanks so much for that info, good to know about the TMS spec, I had assumed wrongly about TMS 2.0 already being implemented. We don't have an immediate need for 2.0 in morecantile but likely that would be of interest eventually. I am still very much learning about the TMS spec in general, so it would be very helpful to get your guidance on what else needs to be done to support these IAU CRSs better in morecantile and also in projects downstream of it.

@AndrewAnnex
Copy link
Owner

@jlaura looks like Vincent merged in the first of the two PRs earlier today, there was a bit more work to do with the 2nd one but I think we're now at a much better spot with it and hope it can be merged in soon. Time to come up with more tasks!

@AndrewAnnex
Copy link
Owner

one thing that could be good is to think about the maxzoom parameter a bit. I could imagine that for low res data sets like MOLA we could tune that to lower numbers but we'd want TMS specs usable for HiRISE images also. Maybe it's best to assume a higher zoom level always for each body and assume clients would do the right thing and not attempt to render tiles much finer than the best resolution available.

@jlaura
Copy link
Contributor Author

jlaura commented Dec 13, 2022

@AndrewAnnex 🎉 That's great news!

I think we need to support a really wide range zoom range and then assume that the client or user will do the correct thing. For example, if we have a TMS MOLA base (resolutions is like 463m or so?) and someone wants to tile a long CTX string (50x lines at 6mpp), I think we want them to be able to use the same TMS endpoints.

Next steps on the TMS:

  • Wait for the opengis.net update to happen.
  • Test deploy the TMSes to the endpoint.
  • Test registering the endpoint to stac browser at stac.astrogeology.usgs.gov/browser-dev. That has Europa, Moon, and Mars example.

Are you thinking of other next steps for TMS or in general? I am not super fluent in the topography encoding, but you seemed pretty excited about that. Can I assist there? I've been starting to look at vector tiles and SLD symbology definitions as well. I have a pygeoapi deploy as an AWS lambda that I am trying to use to learn about the vector side of things. COG and STAC standards / best practices are probably right up there too. For example, we have WKT2 codes, but they will not properly encode out of GDAL without some coaxing. 😆 So, lots of areas!

@AndrewAnnex
Copy link
Owner

@jlaura sounds good regarding the wide zoom ranges, I'll do some thinking to come up with a reasonable limit (possibly 24 or 25; for earth zoom level 24 is 0.009m GSD at the equator assuming an equirectangular projection, so smaller bodies would be even finer scale resolution than that, 24 may be more sensible as that would ensure tiles indexes remain within 24bits if that even matters... some things to think about anyways)

Next steps sound good. Once the morecantile prs are merged in (and maybe a new version is released) the changes in this PR can be merged.

I sent you an email before I saw this and had some thoughts about the topography encoding which for now I think are somewhat independent from this until the opengis update happens. My biggest unknown is what else needs to happen to js frontends assuming we get these parts working. I also haven't thought much about the vector side yet, interested to hear why the wkt2 codes wouldn't work with ogr/gdal...

We may need to move the discussion elsewhere on github, maybe as issues on this repo?

@vincentsarago
Copy link

morecantile 3.2.3 is available on pypi https://pypi.org/project/morecantile/3.2.3/ 🥳
thanks a lot @AndrewAnnex 🙏


for tmsp in crss:
# create the tms object
tms = morecantile.TileMatrixSet.custom(**asdict(tmsp))
_d = json.loads(tms.json(exclude_none=True)) # get to pretty printed json
with open(f'./{tmsp.crs.name}_tms.json', 'w') as dst:
with open(f'./{tmsp.identifier}_tms.json', 'w') as dst:
# dump to json
json.dump(_d, dst, indent=4, ensure_ascii=False)
print(f'wrote {dst.name}')

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So this generate script will create a lot of TMS right? I'm kinda worry now that when importing planetcantile we might try to parse too many JSON files before being able to do any operation.

tms = TileMatrixSet.parse_file(src)

Maybe this could be fixed on morecantile side where we only parse the JSON file when we need it 🤷

Copy link
Contributor Author

@jlaura jlaura Dec 13, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@vincentsarago Correct - a whole bunch of JSON files!

A change to morecantile would be a pretty big change to the logic in morecantile.defaults right? Looks like that is currently globing, loading, and ultimately creating a TileMatrixSet object that is used up stream (e.g., by titiler). Is that correct?

If so, we'd want to make sure that users of Morecantile (and upstream) can see the available TMSs without having to parse all the JSON.

  • We could enforce filename mapping to IDs and just read the filenames (this has issues).
  • We could read just a portion of the JSON (this is hacky and has other issues including a giant WHY).
  • We could read all of them using a library like msgspec which is supposed to be significantly faster. (Caveat, we have lots of file handles to grab for small JSON files, so this would need benchmarking first. Could we go to a single TMS file with an array of TMS objects?).
  • If we could advertise properly, we could read on use and then keep that TMS in memory, so cache.

I guess, I would really want to make sure that we can advertise the availability of the TMS while still maintain good performance? Maybe I'm worried about the wrong thing though!

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I did start a PR to show how it could be developmentseed/morecantile#94

I don't think this is too much of a change. in planetcantile it will look like:

import importlib.resources as importlib_resources
from morecantile.defaults import default_tms, TileMatrixSets

planet_tms_jsons = importlib_resources.files('planetcantile.data').glob('*.json')
default_tms.update(
    {p.stem: p for p in planet_tms_jsons}
)
tms = TileMatrixSets(copy(default_tms))

# Or to only use TMS from planetcantile
import importlib.resources as importlib_resources
from morecantile import TileMatrixSets

planet_tms_jsons = importlib_resources.files('planetcantile.data').glob('*.json')
default_tms = {p.stem: p for p in planet_tms_jsons}
tms = TileMatrixSets(copy(default_tms))

The other solution would be to simply use TILEMATRIXSET_DIRECTORY environment variable to tell morecantile to look for json document within a specific directory. https://github.com/developmentseed/morecantile/blob/main/morecantile/defaults.py#L16-L18

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We could enforce filename mapping to IDs and just read the filenames (this has issues).

Well right now morecantile assume this 🙈

I guess, I would really want to make sure that we can advertise the availability of the TMS while still maintain good performance? Maybe I'm worried about the wrong thing though!

Right now (with the PR) we keep the advertising part (but using the file name).

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well right now morecantile assume this 🙈
😆 This is good! And we can use the environment variable on the titiler side to offer these via our lambda deploy. I might hack the check to opengis.net out of a version that I have deployed (only to be able to test until they update) and then get that URL shared to you both for verification.

Knowing about that env variable and having the code snippet is AWESOME! 🌟

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am worried about projects that use morecantile's default tms set without allowing overrides. Yes it's possible to make PRs against them to accept a tms set as a parameter but my suggestion could be a very cool mechanism to support

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

err "like" it but not exactly it, also not working currently

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am worried about projects that use morecantile's default tms set without allowing overrides.

I'm not sure to understand this

Morecantile is a tool that lets you define and use TMS grid, having a set of defaults in the project is a bonus but isn't the goal of the project.

Yes it's possible to make PRs against them to accept a tms set as a parameter but my suggestion could be a very cool mechanism to support

I'm open for suggestion but right now I have no idea how it could look 🤷

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

err "like" it but not exactly it, also not working currently

what is not working?

@AndrewAnnex
Copy link
Owner

AndrewAnnex commented Dec 13, 2022 via email

@AndrewAnnex
Copy link
Owner

I am going to try to merge this PR and clean up things a bit to ensure that the work is captured. I've saved a PDF of the page saving all of the discussions above just in case but will leave the pr up to allow the conversation to continue a bit here until it needs to move on to other prs/issues

@AndrewAnnex AndrewAnnex merged commit 49e2ab1 into AndrewAnnex:main Dec 13, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants