start of adding projections

[AndrewAnnex]

adds support for projections including equirectangular, north & south polar projections.

includes smarts to compute mostly correct bounds in meters for projections

seeming found issues with ographic projections where the always_xy parameter is not respected when computing bounds. Need to look into that to figure out if it's a bug or unsupported in proj/pyproj

I had been using two saturn CRSS (69910 and 69911) to look into this. which are:

PROJCRS["Saturn (2015) - Sphere / Ocentric/ Equirectangular, clon = 0",
    BASEGEOGCRS["Saturn (2015) - Sphere / Ocentric",
        DATUM["Saturn (2015) - Sphere",
        ELLIPSOID["Saturn (2015) - Sphere", 60268000, 0,
        LENGTHUNIT["metre", 1, ID["EPSG", 9001]]]],
        PRIMEM["Reference Meridian", 0,
            ANGLEUNIT["degree", 0.0174532925199433, ID["EPSG", 9122]]],
        ID["IAU", 69900, 2015]],
    CONVERSION["Equirectangular, clon = 0",
        METHOD["Equidistant Cylindrical",
            ID["EPSG", 1028]],
        PARAMETER["Latitude of 1st standard parallel", 0,
            ANGLEUNIT["degree",0.0174532925199433,ID["EPSG", 9122]],
            ID["EPSG", 8823]],
        PARAMETER["Longitude of natural origin", 0,
            ANGLEUNIT["degree",0.0174532925199433,ID["EPSG", 9122]],
            ID["EPSG", 8802]],
        PARAMETER["False easting", 0,
            LENGTHUNIT["metre",1,ID["EPSG", 9001]],
            ID["EPSG", 8806]],
        PARAMETER["False northing", 0,
            LENGTHUNIT["metre",1,ID["EPSG", 9001]],
            ID["EPSG", 8807]]],
    CS[Cartesian, 2],
        AXIS["Easting (E)", east,
            ORDER[1],
            LENGTHUNIT["metre", 1]],
        AXIS["Northing (N)", north,
            ORDER[2],
            LENGTHUNIT["metre", 1]],
    ID["IAU", 69910, 2015]]

PROJCRS["Saturn (2015) / Ographic/ Equirectangular, clon = 0",
    BASEGEOGCRS["Saturn (2015) / Ographic",
        DATUM["Saturn (2015)",
        ELLIPSOID["Saturn (2015)", 60268000, 10.2079945799458,
        LENGTHUNIT["metre", 1, ID["EPSG", 9001]]]],
        PRIMEM["Reference Meridian", 0,
            ANGLEUNIT["degree", 0.0174532925199433, ID["EPSG", 9122]]],
        ID["IAU", 69901, 2015]],
    CONVERSION["Equirectangular, clon = 0",
        METHOD["Equidistant Cylindrical",
            ID["EPSG", 1028]],
        PARAMETER["Latitude of 1st standard parallel", 0,
            ANGLEUNIT["degree",0.0174532925199433,ID["EPSG", 9122]],
            ID["EPSG", 8823]],
        PARAMETER["Longitude of natural origin", 0,
            ANGLEUNIT["degree",0.0174532925199433,ID["EPSG", 9122]],
            ID["EPSG", 8802]],
        PARAMETER["False easting", 0,
            LENGTHUNIT["metre",1,ID["EPSG", 9001]],
            ID["EPSG", 8806]],
        PARAMETER["False northing", 0,
            LENGTHUNIT["metre",1,ID["EPSG", 9001]],
            ID["EPSG", 8807]]],
    CS[Cartesian, 2],
        AXIS["Westing (W)", west,
            ORDER[1],
            LENGTHUNIT["metre", 1]],
        AXIS["Northing (N)", north,
            ORDER[2],
            LENGTHUNIT["metre", 1]],
    ID["IAU", 69911, 2015]]

The bounds in longitude/latitude should always be either (-180, -90, 180, 90) or (0, -90, 360, 90) depending if clon=0 or not. But when given a sphere vs ellipsoid with a high flattening proj or pyproj doesn't seem to respect the always_xy parameter requiring a swap of x and y in that case.

This could be working as expected, or maybe it's a known issue in proj/pyproj but I don't fully understand the issue yet.

[AndrewAnnex]

@jlaura when you get the chance can you take a look at this and let me know what you think? I think the PR is so big it's breaking the ui preventing me to request reviewers...

[AndrewAnnex]

@thareUSGS, this may also be something you've run into before?

[AndrewAnnex]

@jlaura I wasn't asking you to look at the json files themselves, I was specifically asking about the axis ordering issue and if you had seen it before. It seems like the only thing different about the CRSs is that the 1st axis is a westing rather than an easting. The possibly relevant bit of code in morecantile is https://github.com/developmentseed/morecantile/blob/6334c212e304fd2feec95da34af4fedf2cf994a0/morecantile/models.py#L87-L89 although the crss I posted above seem to both return false correctly. I noticed my note about the clon, maybe I had a typo there...

[AndrewAnnex]

oh I forgot that I was seeing this as a potential issue with pyproj itself there maybe something wrong with transform_bounds...

[AndrewAnnex]

@jlaura okay I think that change resolved the issue fully with the axis ordering, its a little strange to me still that the axis ordering of the geodetic_crs can be so different from the actual wkt text defining the projection though.

I think a good next step I could use from you is a sanity spot-check on some of these. I think the most important details are in the first 40 lines to ensure that the bounding boxes and topLeftCorner positions are correct given the expectations of the WKT.

for example for IAU:49915 I see

{
    "type": "TileMatrixSetType",
    "title": "Mars (2015) - Sphere / Ocentric/ Equirectangular, clon = 180",
    "identifier": "IAU_49915_2015",
    "supportedCRS": "http://www.opengis.net/def/crs/IAU/2015/49915",
    "boundingBox": {
        "type": "BoundingBoxType",
        "crs": "http://www.opengis.net/def/crs/IAU/2015/49915",
        "lowerCorner": [
            -10669445.554195119,
            0.0
        ],
        "upperCorner": [
            10464263.908922138,
            5334722.7770975595
        ]
    },
    "tileMatrix": [
        {
            "type": "TileMatrixType",
            "identifier": "0",
            "scaleDenominator": 147417058.19696748,
            "topLeftCorner": [
                -10669445.554195119,
                5334722.7770975595
            ],
            "tileWidth": 256,
            "tileHeight": 256,
            "matrixWidth": 2,
            "matrixHeight": 1
        },
        {
            "type": "TileMatrixType",
            "identifier": "1",
            "scaleDenominator": 73708529.09848374,
            "topLeftCorner": [
                -10669445.554195119,
                5334722.7770975595
            ],
            "tileWidth": 256,

that lowerCorner 0 I think is a little suspect no?

[AndrewAnnex]

I've made a jupyter notebook to begin looking into this, I think the checks I am doing inside this pr may be doing more harm than good, this may actually be a bug in Pyproj

https://gist.github.com/AndrewAnnex/6b1595f17a892d074902cd33ff833a2d