Playing around with dynamic specifications from wkt2 strings

jlaura commented 1 year ago

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

AndrewAnnex commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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

Validating these through morecantile+pydantic fails because they lack, naturally, EPSG codes. Code below to reproduce. I reached out to the PoC on opengis.net to see about adding the IAU codes. 📞

    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.

```python
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 commented 1 year ago

@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:

Registering the custom TMSs is already a part of planetcantile https://github.com/AndrewAnnex/planetcantile/blob/5ea2577f5dc4a3bc91b0443ef0633a5f89b15e03/planetcantile/defaults.py#L1-L16 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.
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.
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 commented 1 year ago

@jlaura this looks relevant to our discussion https://github.com/opengeospatial/NamingAuthority/issues/212

jlaura commented 1 year ago

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 commented 1 year ago

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 (https://github.com/developmentseed/morecantile/pull/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 commented 1 year ago

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! https://github.com/developmentseed/morecantile/pull/93)

jlaura commented 1 year ago

😆 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 commented 1 year ago

👋 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 commented 1 year ago

@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 commented 1 year ago

@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 commented 1 year ago

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 commented 1 year ago

@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 commented 1 year ago

@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 commented 1 year ago

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

AndrewAnnex commented 1 year ago

Different file in the repository (planetcantile/defaults.py) , no longer broken due to recent pushes I made to fix- Dr. Andrew AnnexOn Dec 13, 2022, at 2:07 PM, Vincent Sarago @.***> wrote: @vincentsarago commented on this pull request.

In planetcantile/data/generate.py:

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}')

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

what is not working?

—Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you were mentioned.Message ID: @.***>

AndrewAnnex commented 1 year ago

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 / planetcantile

Playing around with dynamic specifications from wkt2 strings #2

Option 1: Hack the CRSes in the JSON files.

Option 2: Register a custom TMS to morecantile

create the tms object

dump to json