When using NZTM TileMatrixSet the numbers appear to be off.

import morecantile

nztm = morecantile.tms.get('NZTM2000')
print(nztm.xy_bounds(morecantile.Tile(1, 2, 0)))

returns

CoordsBbox(xmin=12293760.0, ymin=-7881280.0, xmax=14587520.0, ymax=-5587520.0)

When viewing this tile inside of QGIS the top left point is { x: 1293760, y: 5412480 }

I believe this is caused by https://github.com/developmentseed/morecantile/blob/master/morecantile/models.py#L305 as NZTM is defined as [y, x] rather than [x, y]

Hello, since the version bump to 3.1.0 we are getting the following output from pip. It seems like there's an issue with the package? Everything was fine under 3.0.5. Thanks for your help.

#!/bin/bash -eo pipefail
sudo pip3 install cogeo-mosaic
Collecting cogeo-mosaic
  Downloading cogeo-mosaic-4.0.0.tar.gz (23 kB)
  Preparing metadata (setup.py) ... - done
Collecting attrs
  Downloading attrs-21.4.0-py2.py3-none-any.whl (60 kB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 60.6/60.6 KB 32.0 MB/s eta 0:00:00

Collecting cachetools
  Downloading cachetools-5.0.0-py3-none-any.whl (9.1 kB)
Collecting httpx
  Downloading httpx-0.22.0-py3-none-any.whl (84 kB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 84.2/84.2 KB 39.5 MB/s eta 0:00:00

Collecting mercantile
  Downloading mercantile-1.2.1-py3-none-any.whl (14 kB)
Collecting morecantile<4.0,>=3.0
  Downloading morecantile-3.1.0.tar.gz (25 kB)
  Preparing metadata (setup.py) ... - error
  error: subprocess-exited-with-error
  
  × python setup.py egg_info did not run successfully.
  │ exit code: 1
  ╰─> [36 lines of output]
      Traceback (most recent call last):
        File "<string>", line 36, in <module>
        File "<pip-setuptools-caller>", line 34, in <module>
        File "/tmp/pip-install-e4svnxu4/morecantile_2e4a53c9883345079bf39a2a4ca2e893/setup.py", line 45, in <module>
          "morecantile=morecantile.scripts.cli:cli",
        File "/usr/lib/python3/dist-packages/setuptools/__init__.py", line 145, in setup
          return distutils.core.setup(**attrs)
        File "/usr/lib/python3.7/distutils/core.py", line 121, in setup
          dist.parse_config_files()
        File "/usr/lib/python3/dist-packages/setuptools/dist.py", line 705, in parse_config_files
          ignore_option_errors=ignore_option_errors)
        File "/usr/lib/python3/dist-packages/setuptools/config.py", line 120, in parse_configuration
          meta.parse()
        File "/usr/lib/python3/dist-packages/setuptools/config.py", line 425, in parse
          section_parser_method(section_options)
        File "/usr/lib/python3/dist-packages/setuptools/config.py", line 398, in parse_section
          self[name] = value
        File "/usr/lib/python3/dist-packages/setuptools/config.py", line 183, in __setitem__
          value = parser(value)
        File "/usr/lib/python3/dist-packages/setuptools/config.py", line 513, in _parse_version
          version = self._parse_attr(value, self.package_dir)
        File "/usr/lib/python3/dist-packages/setuptools/config.py", line 348, in _parse_attr
          module = import_module(module_name)
        File "/usr/lib/python3.7/importlib/__init__.py", line 127, in import_module
          return _bootstrap._gcd_import(name[level:], package, level)
        File "<frozen importlib._bootstrap>", line 1006, in _gcd_import
        File "<frozen importlib._bootstrap>", line 983, in _find_and_load
        File "<frozen importlib._bootstrap>", line 967, in _find_and_load_unlocked
        File "<frozen importlib._bootstrap>", line 677, in _load_unlocked
        File "<frozen importlib._bootstrap_external>", line 728, in exec_module
        File "<frozen importlib._bootstrap>", line 219, in _call_with_frames_removed
        File "/tmp/pip-install-e4svnxu4/morecantile_2e4a53c9883345079bf39a2a4ca2e893/morecantile/__init__.py", line 14, in <module>
          from .defaults import tms  # noqa
        File "/tmp/pip-install-e4svnxu4/morecantile_2e4a53c9883345079bf39a2a4ca2e893/morecantile/defaults.py", line 8, in <module>
          import attr
      ModuleNotFoundError: No module named 'attr'
      [end of output]
  
  note: This error originates from a subprocess, and is likely not a problem with pip.
error: metadata-generation-failed

× Encountered error while generating package metadata.
╰─> See above for output.

note: This is an issue with the package mentioned above, not pip.
hint: See above for details.


Exited with code exit status 1

Tile matrix sets have the notion of X/Y/Z coordinates so technically we can generate a quadkey. However not all tile matrix sets have 4 children per parent which breaks the parent/child hierarchy. A good example of this is Canadian NAD83 LCC.

This all being said, I do think there are many benefits to supporting quadkeys:

• serialization/deserialization of tiles
• server side caching (generating unique per-tile cache keys)
• fast parent/child lookups

I think this comes down to whether or not morecantile will support tile matrices which don't decimate by a factor of 4 (width or height is 2**zoom). From the README, docstrings, and general usage my thought is that it shouldn't, but that isn't entirely clear.

in #58 we are proposing a major change to this library.

When we started the development of morecantile we choose to use rasterio only to avoid adding requirements like pyproj (because morecantile was mostly made to work with rio-tiler). In addition for a long time rasterio wheels were stuck at proj 6 while pyproj were already using proj>=7, so we wanted to avoid using coordinates coming from 2 different proj version.

Now that rasterio wheels uses proj>=7, if think we are in a better place to switch to pyproj and thus remove rasterio from this libraries.

The main motivation is trying to avoid some weird behaviour with coordinate re-projection which seems to be better handled by proj (with rasterio, you use GDAL which use proj while with pyproj you directly use Proj).

The PR is ready and the CI passes (still trying to match mercantile tests)

What it means for libraries using morecantile

Mostly nothing. The biggest breaking change is the type of CRS being returned by morecantile.TileMatrixSet.crs property, which is now of type pyproj.CRS instead of rasterio.crs.CRS. This is used in rio-tiler for example (https://github.com/cogeotiff/rio-tiler/search?l=Python&q=tms.crs). To help with the transition we are adding morecantile.TileMatrixSet.rasterio_crs property which will return a rasterio.crs.CRS. This method will only be available is rasterio is installed.

import morecantile

# before
tms = morecantile.tms.get("WebMercatorQuad")
crs = tms.crs

# now 
tms = morecantile.tms.get("WebMercatorQuad")
crs = tms.rasterio_crs

If no negative feedback, we'll do a pre-release next week

I'm trying to use this library to create TileMatrixSet definitions (ultimately for cogeo-mosaic, titiler and friends) for non-Earth (Mars) coordinate reference systems. Most aspects of the math are exactly the same, but recent Proj has gotten strict about ellipsoids. This is a good thing for most users, but the "Cannot find coordinate operations from 'EPSG:4326' to XXX" has become a source of dread for planetary scientists ;)

I've been able to successfully adapt the TileMatrixSet definition to mostly work for the "Mars 2000" CRS, but I had to reimplement more than I ideally would to simply change the CRS:

MARS2000_SPHERE = CRS.from_dict(
    {"proj": "longlat", "R": 3396190, "no_defs": True}
)

class MarsTMS(TileMatrixSet):
    _to_mars2000 = PrivateAttr()
    _from_mars2000 = PrivateAttr()

    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)

        # Mars 2000 sphere can be used because lat/lon are defined
        # as planetocentric on Mars
        self._to_mars2000 = Transformer.from_crs(
            self.supportedCRS, MARS2000_SPHERE, always_xy=True
        )
        self._from_mars2000 = Transformer.from_crs(
            MARS2000_SPHERE, self.supportedCRS, always_xy=True
        )

    @property
    def bbox(self):
        """Return TMS bounding box in WGS84."""
        bbox = self._to_mars2000.transform_bounds(*self.xy_bbox, densify_pts=21)
        return BoundingBox(*bbox)

    def xy(self, lng: float, lat: float, truncate=False) -> Coords:
        """Transform longitude and latitude coordinates to TMS CRS."""
        if truncate:
            lng, lat = truncate_lnglat(lng, lat)

        inside = point_in_bbox(Coords(lng, lat), self.bbox)
        if not inside:
            warnings.warn(
                f"Point ({lng}, {lat}) is outside TMS bounds {list(self.bbox)}.",
                PointOutsideTMSBounds,
            )

        x, y = self._from_mars2000.transform(lng, lat)

        return Coords(x, y)

The lengthy redefinition is needed to replace the _to_wgs84 and _from_wgs84 definitions. Ideally, these could be replaced by _to_geographic and _from_geographic methods, with the geographic coordinate system being user-supplied (WGS84 by default). This could also help with other planets and other reasons people might define a custom ellipsoid.

If you think I'm missing a simpler way to achieve this, please let me know. If not, I'm happy to supply a PR.

Dear All,

I haven't heard any feedback since the release of the 3.0.0a0 release so I think we can go ahead and release 3.0.0 early next week.

Regards,

Vincent

in #22 an user wanted to used a custom CRS defined from a proj4 string

import morecantile 
from rasterio.crs import CRS 

crs = CRS.from_proj4("+proj=stere +lat_0=90 +lon_0=0 +k=2 +x_0=0 +y_0=0 +R=3396190 +units=m +no_defs") 
extent = [-13584760.000,-13585240.000,13585240.000,13584760.000] 
tms = morecantile.TileMatrixSet.custom(extent, crs, identifier="MarsNPolek2MOLA5k")  

While the TMS creation works, the internal representation of the CRS won't work because morecantile is trying to create an URI type

tms.supportedCRS                                                                                                                         
>  'http://www.opengis.net/def/crs/EPSG/0/None'

Because

crs = CRS.from_proj4("+proj=stere +lat_0=90 +lon_0=0 +k=2 +x_0=0 +y_0=0 +R=3396190 +units=m +no_defs") 
crs.to_epsg()
> None

I think we should default the tms.crs to the custom one instead of doing inputCRS -> toEPSG -> toURI -> toCRS

cc @cirquelar

I'm trying to create a custom tile matrix based on the example:

import rasterio, morecantile from rasterio.crs import CRS crs = CRS.from_proj4("+proj=stere +lat_0=90 +lon_0=0 +k=2 +x_0=0 +y_0=0 +R=3396190 +units=m +no_defs") extent = [-13584760.000,-13585240.000,13585240.000,13584760.000] tms = morecantile.TileMatrixSet.custom(extent, crs, indentifier="MarsNPolek2MOLA5k") Traceback (most recent call last): File "", line 1, in TypeError: custom() got an unexpected keyword argument 'indentifier'

if I remove 'indentifier':

tms = morecantile.TileMatrixSet.custom(extent, crs, "MarsNPolek2MOLA5k") Traceback (most recent call last): File "", line 1, in File "/Users/fcalef/miniconda3/envs/py37/lib/python3.7/site-packages/morecantile/models.py", line 157, in custom width / (tile_width * matrix_scale[0]) / 2.0 ** zoom, TypeError: unsupported operand type(s) for /: 'float' and 'str'

I don't see any obvious syntax mistakes.

grabs the authority, version and code from the CRS directly to encode the URIs. Related to #92 but can be merged separately. I'm not aware of any issues this would cause but I am open to hearing about them

The referenced specification only considers Meters and Degrees. Would it be possible to consider other units? In the United States were have the State Plane Coordinate System with either defined units in "US Survey Foot" or "International Foot".

US Survey Foot = 0.3048006096 International Foot = 0.3048

I hard-coded the else statement to the US Survey Foot for EPSG:2276 data that I am working with just as a test. And I am very happy to have found success in properly tiling locally referenced imagery!

I don't know Python, but I could take a crack at adding a case statement for the defined units if there is no objection.

https://github.com/developmentseed/morecantile/blob/17020c555d4f09e3027cab755cf0f74744106774/morecantile/utils.py#L55

From note g in http://docs.opengeospatial.org/is/17-083r2/17-083r2.html#table_2: If the CRS uses meters as units of measure for the horizontal dimensions, then metersPerUnit=1; if it has degrees, then metersPerUnit=2pa/360 (a is the Earth maximum radius of the ellipsoid).

The mercantile.tile and morecantile.tms.get("WebMercatorQuad").tile functions are returning slightly different values given the same lat/lon coordinate. This is because mercantile uses slightly more precision when defining the bounds of the Web Mercator tiling scheme. I updated the bounds of the TMS defined in morecantile/data to match the return of mercantile.xy_bounds(mercantile.Tile(0,0,0)). You'll notice that the first build fails, and tests are fixed when making this change.

This change means our WebMercatorQuad deviates slightly from what is defined by OGC (https://docs.opengeospatial.org/is/17-083r2/17-083r2.html#62) but I think its worth it to align with mercantile.

I also had to move the antimeredian test case into its own function, for some reason it fails when its included in the larger test case 🤷

NGA has a "non-web-Mercator" tile grid document ("Map Projections for Tiled Raster Graphics, V1.0, 24 April 2015") available from https://nsgreg.nga.mil/doc/view?i=1810

Its intended for GeoPackage (although application is a bit more general). Possibly useful to incorporate.