diff --git a/.flake8 b/.flake8 index b741c9d..dd7c888 100644 --- a/.flake8 +++ b/.flake8 @@ -1,4 +1,20 @@ [flake8] max-line-length = 100 -ignore = E261,E501,E731,W503 +ignore = + # Ignore line break before binary operator + W503, +per-file-ignores = + __init__.py: + # Ignore unused imports in __init__.py files + F401, + # Ignore 'from module import *' in __init__.py files + F403, + test_*.py: + # Ignore line too long in tests + E501, + # Ignore pydocstyle in tests + D100,D101,D102,D103 + conftest.py: + # Ignore line too long in tests setup + E501, exclude = .ipynb_checkpoints diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 5f495b4..69eaf50 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -15,7 +15,7 @@ jobs: strategy: matrix: os: [macos-latest, ubuntu-latest, windows-latest] - python-version: [3.6, 3.7, 3.8, 3.9] + python-version: [3.6, 3.7, 3.8, 3.9, '3.10'] steps: - uses: actions/checkout@v2 @@ -45,38 +45,3 @@ jobs: run: nox -s pylint - name: Run pydocstyle run: nox -s pydocstyle - - docs: - name: Build and push docs - runs-on: ubuntu-latest - if: github.ref == 'refs/heads/main' - steps: - - uses: actions/checkout@v2 - with: - fetch-depth: 0 # To fetch all refs - - name: Set up Python - uses: actions/setup-python@v1 - with: - python-version: 3.8 - - name: Install nox - run: pip install nox - - name: Build docs - run: nox -s docs -- -d /tmp/sphinx-doctrees - - name: Copy built docs and checkout - run: | - rm -rf __pycache__ **/__pycache__ .doctrees **/.doctrees .buildinfo **/.buildinfo - cp -r docs/build /tmp/sphinx-build - git clean -fx - git checkout gh-pages - git rm -r . # Clean branch to ensure removed files don't persist - - name: Commit and push docs - env: - GIT_WORK_TREE: /tmp/sphinx-build - run: | - git status - git add . - git status - git config --global user.name "${GITHUB_ACTOR}" - git config --global user.email "${GITHUB_ACTOR}@users.noreply.github.com" - git commit --allow-empty --message "Update docs to ${GITHUB_SHA} from ${GITHUB_REF}" - git push diff --git a/.gitignore b/.gitignore index 6aac9ca..f62043d 100644 --- a/.gitignore +++ b/.gitignore @@ -85,4 +85,4 @@ dmypy.json /docs/build # Autogenerated API docs -/docs/source/api +/docs/api diff --git a/.readthedocs.yaml b/.readthedocs.yaml new file mode 100644 index 0000000..32ecbe8 --- /dev/null +++ b/.readthedocs.yaml @@ -0,0 +1,10 @@ +version: 2 +sphinx: + configuration: docs/conf.py +python: + version: 3.7 + install: + - method: pip + path: . + extra_requirements: + - dev diff --git a/README.rst b/README.rst index c6cf5ab..e79db73 100644 --- a/README.rst +++ b/README.rst @@ -17,6 +17,10 @@ library. :alt: CI Status :target: https://github.com/aaliddell/s2cell/actions +.. image:: https://readthedocs.org/projects/s2cell/badge/?version=latest + :alt: Documentation Status + :target: https://docs.s2cell.aliddell.com/en/latest + .. image:: https://img.shields.io/github/license/aaliddell/s2cell :alt: License :target: https://github.com/aaliddell/s2cell diff --git a/docs/source/annotated_source.rst b/docs/annotated_source.rst similarity index 89% rename from docs/source/annotated_source.rst rename to docs/annotated_source.rst index 4609437..45dcb17 100644 --- a/docs/source/annotated_source.rst +++ b/docs/annotated_source.rst @@ -13,6 +13,6 @@ system that are useful for understanding the steps to convert from one represent The full source of the library is included below. -.. literalinclude:: ../../s2cell/__init__.py +.. literalinclude:: ../s2cell/s2cell.py :language: python3 - :lines: 25- + :lines: 22- diff --git a/docs/source/changelog.rst b/docs/changelog.rst similarity index 89% rename from docs/source/changelog.rst rename to docs/changelog.rst index 50281c1..38662d5 100644 --- a/docs/source/changelog.rst +++ b/docs/changelog.rst @@ -5,6 +5,14 @@ Changelog ========= +1.5.0 +----- + +- Moved main source out of ``__init__.py`` +- Moved docs hosting to Read the Docs +- Updated docs to improve readability + + 1.4.0 ----- diff --git a/docs/conf.py b/docs/conf.py index 8f2a27a..317250c 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -1,4 +1,5 @@ import datetime +import os # -- Project information ----------------------------------------------------- @@ -20,20 +21,20 @@ extensions = [ # Internal 'sphinx.ext.autodoc', - 'sphinx.ext.githubpages', 'sphinx.ext.napoleon', # External + 'notfound.extension', 'sphinx_sitemap', ] # Add any paths that contain templates here, relative to this directory. -templates_path = [] +templates_path = ['templates'] # List of patterns, relative to source directory, that match files and # directories to ignore when looking for source files. # This pattern also affects html_static_path and html_extra_path. -exclude_patterns = ['.ipynb_checkpoints', '**/.ipynb_checkpoints'] +exclude_patterns = ['.ipynb_checkpoints', '**/.ipynb_checkpoints', 'api/s2cell.rst'] # Code highlighting pygments_style = 'monokai' @@ -41,59 +42,88 @@ # Enable figure numbering numfig = True +# Setup root doc +root_doc = 'index' + + +# -- Run sphinx-apidoc ------------------------------------------------------- + +def run_apidoc(_): + import sphinx.ext.apidoc + + docs_path = os.path.dirname(__file__) + apidoc_path = os.path.join(docs_path, 'api') + module_path = os.path.join(docs_path, '..', 's2cell') + + sphinx.ext.apidoc.main([ + '--no-toc', + '--force', + '--separate', + '-o', apidoc_path, + module_path + ]) + + +def setup(app): + app.connect('builder-inited', run_apidoc) + # -- Options for HTML output ------------------------------------------------- # The theme to use for HTML and HTML Help pages. See the documentation for # a list of builtin themes. # -import sphinx_redactor_theme -html_theme = 'sphinx_redactor_theme' -html_theme_path = [sphinx_redactor_theme.get_html_theme_path()] +html_theme = 'furo' + +# Title +html_title = 's2cell' # Logo and favicon -html_logo = 'source/_static/logo.min.svg' -html_favicon = 'source/_static/logo-64.png' +html_logo = 'static/logo.min.svg' +html_favicon = 'static/logo-64.png' # Base URL for docs -# Used to generate CNAME file html_baseurl = 'https://docs.s2cell.aliddell.com' # Extra vars to provide to templating html_context = { 'baseurl': html_baseurl, - 'icon_png': 'logo.png' + 'icon_png': 'logo.png', + 'absolute_icon_png': html_baseurl + '/en/latest/_static/logo.png' # Used by meta tags } # Add any paths that contain custom static files (such as style sheets) here, # relative to this directory. They are copied after the builtin static files, # so a file named "default.css" will overwrite the builtin "default.css". -html_static_path = ['source/_static'] +html_static_path = ['static'] # Theme options html_theme_options = { - 'display_version': True, + 'sidebar_hide_name': True, + 'light_css_variables': { + 'color-brand-primary': '#de6b00', + 'color-brand-content': '#de6b00', + }, + 'dark_css_variables': { + 'color-brand-primary': '#de6b00', + 'color-brand-content': '#de6b00', + }, } # Disable footer html_show_sphinx = False -# Template overrides -templates_path = ['_templates'] - # Add CSS files -html_css_files = [ - 'pygments.css', # Manually add pygments.css, since sphinx_redactor_theme does not - 'customise.css', -] +html_css_files = [] # Extra files to include html_extra_path = [ - 'source/robots.txt', + 'robots.txt', ] # Sitemap options -sitemap_url_scheme = '{lang}{link}' +sitemap_filename = "sitemap-override.xml" # RTD generates sitemap with not much in it... +sitemap_url_scheme = "{link}" # -- Options for autodoc ----------------------------------------------------- diff --git a/docs/source/index.rst b/docs/index.rst similarity index 67% rename from docs/source/index.rst rename to docs/index.rst index 16ec148..881f79c 100644 --- a/docs/source/index.rst +++ b/docs/index.rst @@ -2,15 +2,15 @@ :description: Minimal Python S2 cell ID, S2 token and lat/lon conversion library :keywords: S2, S2 Geometry, s2cell, S2 cell, cell ID, S2 token, token, Python -.. include:: ../../README.rst +.. include:: ../README.rst .. toctree:: :hidden: Overview - s2_concepts - API Reference - useful_s2_links + Concepts + API Reference + Useful Links annotated_source changelog diff --git a/docs/source/robots.txt b/docs/robots.txt similarity index 100% rename from docs/source/robots.txt rename to docs/robots.txt diff --git a/docs/source/s2_concepts.rst b/docs/s2_concepts.rst similarity index 99% rename from docs/source/s2_concepts.rst rename to docs/s2_concepts.rst index ae16f08..d700772 100644 --- a/docs/source/s2_concepts.rst +++ b/docs/s2_concepts.rst @@ -34,7 +34,7 @@ position generally become close in 2D space). .. _fig_hilbert: -.. figure:: _static/hilbert.svg +.. figure:: static/hilbert.svg :height: 200px :alt: Hilbert Curve Steps 1, 2 and 3 (Qef, Public domain, via Wikimedia Commons) :align: center @@ -144,7 +144,7 @@ the base Hilbert Curve. .. _fig_cube_unwrapped: -.. figure:: _static/cube_unwrapped.svg +.. figure:: static/cube_unwrapped.svg :alt: Cube face mapping :align: center @@ -170,7 +170,7 @@ cube. .. _fig_uv_face_0: -.. figure:: _static/uv_face_0.svg +.. figure:: static/uv_face_0.svg :alt: Face 0 in UV coordinates :align: center @@ -207,7 +207,7 @@ projection. :numref:`fig_st_face_0` shows the ST mapping of the same region as s .. _fig_uv_to_st_projections: -.. figure:: _static/uv_to_st_projections.svg +.. figure:: static/uv_to_st_projections.svg :alt: The three UV to ST projections provided in the reference S2 implementation :align: center @@ -215,7 +215,7 @@ projection. :numref:`fig_st_face_0` shows the ST mapping of the same region as s .. _fig_st_face_0: -.. figure:: _static/st_face_0.svg +.. figure:: static/st_face_0.svg :alt: Face 0 in ST coordinates :align: center diff --git a/docs/source/_static/customise.css b/docs/source/_static/customise.css deleted file mode 100644 index 1a43d82..0000000 --- a/docs/source/_static/customise.css +++ /dev/null @@ -1,166 +0,0 @@ -/* Fix line height of headings */ -h1, h2, h3, h4, h5, .site-content .topic-title, .toctree-wrapper .caption, .auto-toc .caption { - overflow: auto; - line-height: normal; - background: none; -} - -/* Reduce horizontal padding in left sidebar */ -section.site-sidebar { - padding-left: 4%; - padding-right: 4%; -} - -/* Limit max width of main text */ -main.site-main > * { - max-width: 70em; - margin-left: auto; - margin-right: auto; -} - -/* Centre align logo and title in sidebar and reduce its size when too big */ -a.branding-link { - color: black; - font-size: 200%; - font-weight: bold; - padding: 0; - text-align: center; -} - -img.logo { - max-width: 100%; -} - -/* Increase font size in sidebar */ -section.site-nav { - font-size: 100%; -} - -/* Colour hamburger menu on mobile */ -a.navigation-toggle { - color: #de6b00; -} - -/* Colour search box outline */ -.site-searchform [type="text"]:focus { - box-shadow: none; - border-color: #de6b00; -} - -/* Colour links */ -a, .site-content a .literal { - color: #de6b00; -} - -a:hover, .site-content a:hover .literal { - color: unset; -} - -/* Colour code blocks to work with monokai pygments theme */ -div.highlight-bash, div.highlight-default, div.highlight-python3, div.highlight { - background-color: #333; -} - -/* Colour table headings */ -th { - background-color: #de6b00; -} - -/* Colour footer background and text */ -footer.site-footer { - background-color: #de6b00; -} - -footer.site-footer > div > div > p { - color: white; -} - -/* Limit image widths */ -img { - max-width: 100%; -} - -/* Hide prev/next buttons */ -.btn--primary { - display: none; -} - -/* Hide unusable scroll bars on code blocks */ -div.highlight { - overflow-y: hidden; -} - -/* Restyle API reference sections */ -.site-content dl.function { - margin-top: 40px; -} - -.site-content dl:not(.docutils) dt { - border-radius: 0; - background: none; - padding: 0 20px; - margin: 10px 0; -} - -/* Remove colour overloads */ -.branding-link__version, .site-breadcrumbs, .site-breadcrumbs__leaf, .site-sidebar a, -.site-content .literal { - color: unset; -} - -/* Fix how mobile menu behaves to hide element and stop Google complaining */ -@media only screen and (max-width: 980px) { - .site-main { - transition: none; - } - - .site-sidebar { - display: none; - transition: none; - } - - .mobile-menu-open .site-sidebar { - display: block; - } -} - -/* Dark theme */ -@media (prefers-color-scheme: dark) { - body { - background-color: #333; - color: white; - } - - .site-sidebar { - background-color: #333; - } - - .site-main { - background-color: #222; - } - - .navigation-toggle, .mobile-menu-open .navigation-toggle { - background: rgba(0, 0, 0, 0.3); - } - - .site-content .literal { - background-color: #333; - } - - table { - background-color: #444; - } - - table .row-odd td { - background: #333; - } - - td { - border-bottom: 1px solid #222; - } - - .site-searchform [type="text"] { - background-color: #222; - color: unset; - } -} diff --git a/docs/source/_static/cube_unwrapped.png b/docs/static/cube_unwrapped.png similarity index 100% rename from docs/source/_static/cube_unwrapped.png rename to docs/static/cube_unwrapped.png diff --git a/docs/source/_static/cube_unwrapped.svg b/docs/static/cube_unwrapped.svg similarity index 100% rename from docs/source/_static/cube_unwrapped.svg rename to docs/static/cube_unwrapped.svg diff --git a/docs/source/_static/hilbert.svg b/docs/static/hilbert.svg similarity index 100% rename from docs/source/_static/hilbert.svg rename to docs/static/hilbert.svg diff --git a/docs/source/_static/logo-200.png b/docs/static/logo-200.png similarity index 100% rename from docs/source/_static/logo-200.png rename to docs/static/logo-200.png diff --git a/docs/source/_static/logo-64.png b/docs/static/logo-64.png similarity index 100% rename from docs/source/_static/logo-64.png rename to docs/static/logo-64.png diff --git a/docs/source/_static/logo.min.svg b/docs/static/logo.min.svg similarity index 100% rename from docs/source/_static/logo.min.svg rename to docs/static/logo.min.svg diff --git a/docs/source/_static/logo.png b/docs/static/logo.png similarity index 100% rename from docs/source/_static/logo.png rename to docs/static/logo.png diff --git a/docs/source/_static/logo.svg b/docs/static/logo.svg similarity index 100% rename from docs/source/_static/logo.svg rename to docs/static/logo.svg diff --git a/docs/source/_static/st_face_0.png b/docs/static/st_face_0.png similarity index 100% rename from docs/source/_static/st_face_0.png rename to docs/static/st_face_0.png diff --git a/docs/source/_static/st_face_0.svg b/docs/static/st_face_0.svg similarity index 100% rename from docs/source/_static/st_face_0.svg rename to docs/static/st_face_0.svg diff --git a/docs/source/_static/uv_face_0.png b/docs/static/uv_face_0.png similarity index 100% rename from docs/source/_static/uv_face_0.png rename to docs/static/uv_face_0.png diff --git a/docs/source/_static/uv_face_0.svg b/docs/static/uv_face_0.svg similarity index 100% rename from docs/source/_static/uv_face_0.svg rename to docs/static/uv_face_0.svg diff --git a/docs/source/_static/uv_to_st_projections.png b/docs/static/uv_to_st_projections.png similarity index 100% rename from docs/source/_static/uv_to_st_projections.png rename to docs/static/uv_to_st_projections.png diff --git a/docs/source/_static/uv_to_st_projections.svg b/docs/static/uv_to_st_projections.svg similarity index 100% rename from docs/source/_static/uv_to_st_projections.svg rename to docs/static/uv_to_st_projections.svg diff --git a/docs/_templates/layout.html b/docs/templates/base.html similarity index 77% rename from docs/_templates/layout.html rename to docs/templates/base.html index 1acb3c3..57d0523 100644 --- a/docs/_templates/layout.html +++ b/docs/templates/base.html @@ -1,11 +1,11 @@ -{% extends "!layout.html" %} +{% extends "!base.html" %} {% block extrahead %} {%- if meta %} {%- for tag in ('description', 'author', 'keywords') %} {%- if tag in meta %} - + {%- endif %} {%- endfor %} {%- endif %} @@ -15,7 +15,7 @@ {%- endif %} {%- if icon_png %} - + {%- endif %} @@ -25,15 +25,15 @@ {%- if baseurl %} - + {%- endif %} - + {%- if baseurl and icon_png %} - + {%- endif %} - + {%- if meta and 'description' in meta %} - + {%- endif %} diff --git a/docs/source/useful_s2_links.rst b/docs/useful_s2_links.rst similarity index 91% rename from docs/source/useful_s2_links.rst rename to docs/useful_s2_links.rst index ed8dd1b..88b6dc8 100644 --- a/docs/source/useful_s2_links.rst +++ b/docs/useful_s2_links.rst @@ -11,14 +11,15 @@ Your micro 'Awesome S2 Geometry' list. If you have another S2 related link that please `open an Issue `__ or PR. -Concepts --------- +Documentation +------------- -Core concepts of S2 and the S2 cell system. +Core documentation of S2 and the S2 cell system. - `S2 Geometry `__: The S2 Geometry homepage. - `S2 Cells `__: Reference S2 documentation on the S2 cell system. +- `S2 Geometry Concepts `__: The concepts page of this documentation. - `Earth Cube `__: Description of the face cell mapping in the S2 cell system. - `S2 Cell Statistics `__: Details on the sizes @@ -69,6 +70,8 @@ Users of S2 in general, not just of this library. Used by the ``geo_point_to_s2cell`` function. - `BigQuery `__: Used to implement the BigQuery geography functions. +- `ClickHouse `__: Used to implement + geographical data and indexing. - `CockroachDB `__: Used to implement the spatial indexing. - `InfluxDB `__: diff --git a/noxfile.py b/noxfile.py index 1b75015..41e9125 100644 --- a/noxfile.py +++ b/noxfile.py @@ -42,9 +42,6 @@ def build_wheel(session): @nox.session() def docs(session): session.install('.[dev]') - session.run( - 'sphinx-apidoc', '--no-toc', '--force', '--separate', '-o', 'docs/source/api', 's2cell' - ) session.run( 'python', '-m', 'sphinx', '-c', 'docs/', @@ -54,7 +51,7 @@ def docs(session): '-W', # Treat warnings as errors '--keep-going', # When using -W, only exit after all warnings shown #'-j', 'auto', # Generate in parallel - 'docs/source', 'docs/build', + 'docs', 'docs/build', ) @nox.session() diff --git a/pylintrc b/pylintrc index ce708a1..c27656e 100644 --- a/pylintrc +++ b/pylintrc @@ -22,7 +22,7 @@ ignore-patterns= # Use multiple processes to speed up Pylint. Specifying 0 will auto-detect the # number of processors available to use. -jobs=1 +jobs=0 # Control the amount of potential inferred values when inferring a single # object. This can help the performance when dealing with large functions or @@ -60,85 +60,7 @@ confidence= # --enable=similarities". If you want to run only the classes checker, but have # no Warning level messages displayed, use "--disable=all --enable=classes # --disable=W". -disable=print-statement, - parameter-unpacking, - unpacking-in-except, - old-raise-syntax, - backtick, - long-suffix, - old-ne-operator, - old-octal-literal, - import-star-module-level, - non-ascii-bytes-literal, - raw-checker-failed, - bad-inline-option, - locally-disabled, - file-ignored, - suppressed-message, - useless-suppression, - deprecated-pragma, - use-symbolic-message-instead, - apply-builtin, - basestring-builtin, - buffer-builtin, - cmp-builtin, - coerce-builtin, - execfile-builtin, - file-builtin, - long-builtin, - raw_input-builtin, - reduce-builtin, - standarderror-builtin, - unicode-builtin, - xrange-builtin, - coerce-method, - delslice-method, - getslice-method, - setslice-method, - no-absolute-import, - old-division, - dict-iter-method, - dict-view-method, - next-method-called, - metaclass-assignment, - indexing-exception, - raising-string, - reload-builtin, - oct-method, - hex-method, - nonzero-method, - cmp-method, - input-builtin, - round-builtin, - intern-builtin, - unichr-builtin, - map-builtin-not-iterating, - zip-builtin-not-iterating, - range-builtin-not-iterating, - filter-builtin-not-iterating, - using-cmp-argument, - eq-without-hash, - div-method, - idiv-method, - rdiv-method, - exception-message-attribute, - invalid-str-codec, - sys-max-int, - bad-python3-import, - deprecated-string-function, - deprecated-str-translate-call, - deprecated-itertools-function, - deprecated-types-field, - next-method-defined, - dict-items-not-iterating, - dict-keys-not-iterating, - dict-values-not-iterating, - deprecated-operator-function, - deprecated-urllib-function, - xreadlines-attribute, - deprecated-sys-function, - exception-escape, - comprehension-escape +disable=consider-using-f-string # Enable the message, report, category or checker with the given id(s). You can # either give multiple identifier separated by comma (,) or put this option @@ -249,7 +171,7 @@ logging-modules=logging # This flag controls whether inconsistent-quotes generates a warning when the # character used as a quote delimiter is used inconsistently within a module. -check-quote-consistency=no +check-quote-consistency=yes # This flag controls whether the implicit-str-concat should generate a warning # on implicit string concatenation in sequences defined over several lines. diff --git a/s2cell/__init__.py b/s2cell/__init__.py index 1a7faea..22df2a2 100644 --- a/s2cell/__init__.py +++ b/s2cell/__init__.py @@ -14,936 +14,6 @@ """Minimal Python S2 Geometry cell ID, token and lat/lon conversion library.""" -import math -import re -from typing import Optional, Tuple +from .s2cell import * -__version__ = '1.4.0' - - -# -# s2cell exceptions -# - -class InvalidCellID(Exception): - """Exception type for invalid cell IDs.""" - - -class InvalidToken(Exception): - """Exception type for invalid tokens.""" - - -# -# S2 base constants needed for cell mapping -# - -# The maximum level supported within an S2 cell ID. Each level is represented by two bits in the -# final cell ID -_S2_MAX_LEVEL = 30 - -# The maximum value within the I and J bits of an S2 cell ID -_S2_MAX_SIZE = 1 << _S2_MAX_LEVEL - -# The number of bits in a S2 cell ID used for specifying the base face -_S2_FACE_BITS = 3 - -# The number of bits in a S2 cell ID used for specifying the position along the Hilbert curve -_S2_POS_BITS = 2 * _S2_MAX_LEVEL + 1 - -# The maximum value of the Si/Ti integers used when mapping from IJ to ST. This is twice the max -# value of I and J, since Si/Ti allow referencing both the center and edge of a leaf cell -_S2_MAX_SI_TI = 1 << (_S2_MAX_LEVEL + 1) - -# Mask that specifies the swap orientation bit for the Hilbert curve -_S2_SWAP_MASK = 1 - -# Mask that specifies the invert orientation bit for the Hilbert curve -_S2_INVERT_MASK = 2 - -# The number of bits per I and J in the lookup tables -_S2_LOOKUP_BITS = 4 - -# Lookup table for mapping 10 bits of IJ + orientation to 10 bits of Hilbert curve position + -# orientation. Populated later by _s2_init_lookups -_S2_LOOKUP_POS = None - -# Lookup table for mapping 10 bits of Hilbert curve position + orientation to 10 bits of IJ + -# orientation. Populated later by _s2_init_lookups -_S2_LOOKUP_IJ = None - -# Lookup table of two bits of IJ from two bits of curve position, based also on the current curve -# orientation from the swap and invert bits -_S2_POS_TO_IJ = [ - [0, 1, 3, 2], # 0: Normal order, no swap or invert - [0, 2, 3, 1], # 1: Swap bit set, swap I and J bits - [3, 2, 0, 1], # 2: Invert bit set, invert bits - [3, 1, 0, 2], # 3: Swap and invert bits set -] - -# Lookup for the orientation update mask of one of the four sub-cells within a higher level cell. -# This mask is XOR'ed with the current orientation to get the sub-cell orientation. -_S2_POS_TO_ORIENTATION_MASK = [_S2_SWAP_MASK, 0, 0, _S2_SWAP_MASK | _S2_INVERT_MASK] - - -# -# S2 helper functions -# - -def _s2_uv_to_st(component: float) -> float: - """ - Convert S2 UV to ST. - - This is done using the quadratic projection that is used by default for S2. The C++ and Java S2 - libraries use a different definition of the ST cell-space, but the end result in IJ is the same. - The below uses the C++ ST definition. - - See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2coords.h#L317-L320 - - """ - if component >= 0.0: - return 0.5 * math.sqrt(1.0 + 3.0 * component) - return 1.0 - 0.5 * math.sqrt(1.0 - 3.0 * component) - - -def _s2_st_to_uv(component: float) -> float: - """ - Convert S2 ST to UV. - - This is done using the quadratic projection that is used by default for S2. The C++ and Java S2 - libraries use a different definition of the ST cell-space, but the end result in IJ is the same. - The below uses the C++ ST definition. - - See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2coords.h#L312-L315 - - """ - if component >= 0.5: - return (1.0 / 3.0) * (4.0 * component ** 2 - 1.0) - return (1.0 / 3.0) * (1.0 - 4.0 * (1.0 - component) ** 2) - - -def _s2_st_to_ij(component: float) -> int: - """ - Convert S2 ST to IJ. - - The mapping here differs between C++ and Java versions, but the combination of - _st_to_ij(_uv_to_st(val)) is the same for both. The below uses the C++ ST definition. - - See s2geometry/blob/2c02e21040e0b82aa5719e96033d02b8ce7c0eff/src/s2/s2coords.h#L333-L336 - - """ - # The reference implementation does round(_S2_MAX_SIZE * component - 0.5), which is equivalent - # to math.floor(_S2_MAX_SIZE * component) - return int(max(0, min(_S2_MAX_SIZE - 1, math.floor(_S2_MAX_SIZE * component)))) - - -def _s2_si_ti_to_st(component: int) -> float: - """ - Convert S2 Si/Ti to ST. - - This converts an integer in range 0 to _S2_MAX_SI_TI into a float in range 0.0 to 1.0. - - See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2coords.h#L338-L341 - - """ - return (1.0 / _S2_MAX_SI_TI) * component - - -def _s2_face_uv_to_xyz( # pylint: disable=invalid-name - face: int, uv: Tuple[float, float] -) -> Tuple[float, float, float]: - """ - Convert face + UV to S2Point XYZ. - - See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2coords.h#L348-L357 - - Args: - face: The S2 face for the input point. - uv: The S2 face UV coordinates. - - Returns: - The unnormalised S2Point XYZ. - - Raises: - ValueError: If the face is not valid in range 0-5. - - """ - # Face -> XYZ components -> indices with negation: - # 0 -> ( 1, u, v) -> ( /, 0, 1) - # 1 -> (-u, 1, v) -> (-0, /, 1) - # 2 -> (-u, -v, 1) -> (-0, -1, /) - # 3 -> (-1, -v, -u) -> (-/, -1, -0) <- -1 here means -1 times the value in index 1, - # 4 -> ( v, -1, -u) -> ( 1, -/, -0) not index -1 - # 5 -> ( v, u, -1) -> ( 1, 0, -/) - if face == 0: - s2_point = (1, uv[0], uv[1]) - elif face == 1: - s2_point = (-uv[0], 1, uv[1]) - elif face == 2: - s2_point = (-uv[0], -uv[1], 1) - elif face == 3: - s2_point = (-1, -uv[1], -uv[0]) - elif face == 4: - s2_point = (uv[1], -1, -uv[0]) - elif face == 5: - s2_point = (uv[1], uv[0], -1) - else: - raise ValueError('Cannot convert UV to XYZ with invalid face: {}'.format(face)) - - return s2_point - - -def _s2_init_lookups() -> None: - """ - Initialise the S2 lookups in global vars _S2_LOOKUP_POS and _S2_LOOKUP_IJ. - - This generates 4 variations of a 4 level deep Hilbert curve, one for each swap/invert bit - combination. This allows mapping between 8 bits (+2 orientation) of Hilbert curve position and 8 - bits (+2 orientation) of I and J, and vice versa. The new orientation bits read from the mapping - tell us the base orientation of the curve segments within the next deeper level of sub-cells. - - This implementation differs in structure from the reference implementation, since it is - iterative rather than recursive. The end result is the same lookup table. - - See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.cc#L75-L109 - - """ - global _S2_LOOKUP_POS, _S2_LOOKUP_IJ # pylint: disable=global-statement - if _S2_LOOKUP_POS is None or _S2_LOOKUP_IJ is None: # pragma: no branch - # Initialise empty lookup tables - lookup_length = 1 << int(2 * _S2_LOOKUP_BITS + 2) # = 1024 - _S2_LOOKUP_POS = [0] * lookup_length - _S2_LOOKUP_IJ = [0] * lookup_length - - # Generate lookups for each of the base orientations given by the swap and invert bits - for base_orientation in [ - 0, _S2_SWAP_MASK, _S2_INVERT_MASK, _S2_SWAP_MASK | _S2_INVERT_MASK # 0-3 effectively - ]: - # Walk the 256 possible positions within a level 4 curve. This implementation is not - # the fastest since it does not reuse the common ancestor of neighbouring positions, but - # is simpler to read - for pos in range(4 ** 4): # 4 levels of sub-divisions - ij = 0 # Has pattern iiiijjjj, not ijijijij # pylint: disable=invalid-name - orientation = base_orientation - - # Walk the pairs of bits of pos, from most significant to least, getting IJ and - # orientation as we go - for bit_pair_offset in range(4): - # Bit pair is effectively the sub-cell index - bit_pair = (pos >> ((3 - bit_pair_offset) * 2)) & 0b11 - - # Get the I and J for the sub-cell index. These need to be spread into iiiijjjj - # by inserting as bit positions 4 and 0 - ij_bits = _S2_POS_TO_IJ[orientation][bit_pair] - ij = ( # pylint: disable=invalid-name - (ij << 1) # Free up position 4 and 0 from old IJ - | ((ij_bits & 2) << 3) # I bit in position 4 - | (ij_bits & 1) # J bit in position 0 - ) - - # Update the orientation with the new sub-cell orientation - orientation = orientation ^ _S2_POS_TO_ORIENTATION_MASK[bit_pair] - - # Shift IJ and position to allow orientation bits in LSBs of lookup - ij <<= 2 # pylint: disable=invalid-name - pos <<= 2 - - # Write lookups - _S2_LOOKUP_POS[ij | base_orientation] = pos | orientation - _S2_LOOKUP_IJ[pos | base_orientation] = ij | orientation - - -# -# Cell ID <-> Token translation functions -# - -def cell_id_to_token(cell_id: int) -> str: - """ - Convert S2 cell ID to a S2 token. - - Converts the S2 cell ID to hex and strips any trailing zeros. The 0 cell ID token is represented - as 'X' to prevent it being an empty string. - - See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.cc#L204-L220 - - Args: - cell_id: The S2 cell ID integer. - - Returns: - The S2 token string for the S2 cell ID. - - Raises: - TypeError: If the cell_id is not int. - - """ - # Check type - if not isinstance(cell_id, int): - raise TypeError('Cannot convert S2 cell ID from type: {}'.format(type(cell_id))) - - # The zero token is encoded as 'X' rather than as a zero-length string - if cell_id == 0: - return 'X' - - # Convert cell ID to 16 character hex string and strip any implicit trailing zeros - return '{:016x}'.format(cell_id).rstrip('0') - - -def token_to_cell_id(token: str) -> int: - """ - Convert S2 token to S2 cell ID. - - Restores the stripped 0 characters from the token and converts the hex string to integer. - - See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.cc#L222-L239 - - Args: - token: The S2 token string. Can be upper or lower case hex string. - - Returns: - The S2 cell ID for the S2 token. - - Raises: - TypeError: If the token is not str. - InvalidToken: If the token length is over 16. - - """ - # Check input - if not isinstance(token, str): - raise TypeError('Cannot convert S2 token from type: {}'.format(type(token))) - - if len(token) > 16: - raise InvalidToken('Cannot convert S2 token with length > 16 characters') - - # Check for the zero cell ID represented by the character 'x' or 'X' rather than as the empty - # string - if token in ('x', 'X'): - return 0 - - # Add stripped implicit zeros to create the full 16 character hex string - token = token + ('0' * (16 - len(token))) - - # Convert to cell ID by converting hex to int - return int(token, 16) - - -# -# Encode functions -# - -def lat_lon_to_cell_id( # pylint: disable=too-many-locals - lat: float, lon: float, level: int = 30 -) -> int: - """ - Convert lat/lon to a S2 cell ID. - - It is expected that the lat/lon provided are normalised, with latitude in the range -90 to 90. - - Args: - lat: The latitude to convert, in degrees. - lon: The longitude to convert, in degrees. - level: The level of the cell ID to generate, from 0 up to 30. - - Returns: - The S2 cell ID for the lat/lon location. - - Raises: - ValueError: When level is not an integer, is < 0 or is > 30. - - """ - if not isinstance(level, int) or level < 0 or level > _S2_MAX_LEVEL: - raise ValueError('S2 level must be integer >= 0 and <= 30') - - # Populate _S2_LOOKUP_POS on first run. - # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.cc#L75-L109 - # - # This table takes 10 bits of I and J and orientation and returns 10 bits of curve position and - # new orientation - if _S2_LOOKUP_POS is None: # pragma: no cover - _s2_init_lookups() - - # Reuse constant expressions - lat_rad = math.radians(lat) - lon_rad = math.radians(lon) - sin_lat_rad = math.sin(lat_rad) - cos_lat_rad = math.cos(lat_rad) - sin_lon_rad = math.sin(lon_rad) - cos_lon_rad = math.cos(lon_rad) - - # Convert to S2Point - # This is effectively the unit non-geodetic ECEF vector - s2_point = ( - cos_lat_rad * cos_lon_rad, # X - cos_lat_rad * sin_lon_rad, # Y - sin_lat_rad # Z - ) - - # Get cube face - # See s2geometry/blob/2c02e21040e0b82aa5719e96033d02b8ce7c0eff/src/s2/s2coords.h#L380-L384 - # - # The face is determined by the largest XYZ component of the S2Point vector. When the component - # is negative, the second set of three faces is used. - # Largest component -> face: - # +x -> 0 - # +y -> 1 - # +z -> 2 - # -x -> 3 - # -y -> 4 - # -z -> 5 - face = max(enumerate(s2_point), key=lambda p: abs(p[1]))[0] # Largest absolute component - if s2_point[face] < 0.0: - face += 3 - - # Convert face + XYZ to cube-space face + UV - # See s2geometry/blob/2c02e21040e0b82aa5719e96033d02b8ce7c0eff/src/s2/s2coords.h#L366-L372 - # - # The faces are oriented to ensure continuity of curve. - # Face -> UV components -> indices with negation (without divisor, which is always the remaining - # component (index: face % 3)): - # 0 -> ( y, z) -> ( 1, 2) - # 1 -> (-x, z) -> (-0, 2) - # 2 -> (-x, -y) -> (-0, -1) <- -1 here means -1 times the value in index 1, not index -1 - # 3 -> ( z, y) -> ( 2, 1) - # 4 -> ( z, -x) -> ( 2, -0) - # 5 -> (-y, -x) -> (-1, -0) - # - # For a compiled language, a switch statement on face is preferable as it will be more easily - # optimised as a jump table etc; but in Python the indexing method is more concise. - # - # The index selection can be reduced to some bit magic: - # U: 1 - ((face + 1) >> 1) - # V: 2 - (face >> 1) - # - # The negation of the the two components is then selected: - # U: (face in [1, 2, 5]) ? -1: 1 - # V: (face in [2, 4, 5])) ? -1: 1 - uv = ( # pylint: disable=invalid-name - s2_point[1 - ((face + 1) >> 1)] / s2_point[face % 3], # U - s2_point[2 - (face >> 1)] / s2_point[face % 3] # V - ) - if face in (1, 2, 5): - uv = (-uv[0], uv[1]) # Negate U # pylint: disable=invalid-name - if face in (2, 4, 5): - uv = (uv[0], -uv[1]) # Negate V # pylint: disable=invalid-name - - # Project cube-space UV to cell-space ST - # See s2geometry/blob/2c02e21040e0b82aa5719e96033d02b8ce7c0eff/src/s2/s2coords.h#L317-L320 - st = (_s2_uv_to_st(uv[0]), _s2_uv_to_st(uv[1])) # pylint: disable=invalid-name - - # Convert ST to IJ integers - # See s2geometry/blob/2c02e21040e0b82aa5719e96033d02b8ce7c0eff/src/s2/s2coords.h#L333-L336 - ij = (_s2_st_to_ij(st[0]), _s2_st_to_ij(st[1])) # pylint: disable=invalid-name - - # Convert face + IJ to cell ID - # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.cc#L256-L298 - # - # This is done by looking up 8 bits of I and J (4 each) at a time in the lookup table, along - # with two bits of orientation (swap (1) and invert (2)). This gives back 8 bits of position - # along the curve and two new orientation bits for the curve within the sub-cells in the next - # step. - # - # The swap bit swaps I and J with each other - # The invert bit inverts the bits of I and J, which means axes are negated - # - # Compared to the standard versions, we check the required number of steps we need to do for the - # requested level and don't perform steps that will be completely overwritten in the truncation - # below, rather than always doing every step. Each step does 4 bits each of I and J, which is 4 - # levels, so the required number of steps is ceil((level + 2) / 4), when level is > 0. The - # additional 2 levels added are required to account for the top 3 bits (4 before right shift) - # that are occupied by the face bits. - bits = face & _S2_SWAP_MASK # iiiijjjjoo. Initially set by by face - cell_id = face << (_S2_POS_BITS - 1) # Insert face at second most signficant bits - lookup_mask = (1 << int(_S2_LOOKUP_BITS)) - 1 # Mask of 4 one bits: 0b1111 - required_steps = math.ceil((level + 2) / 4) if level > 0 else 0 - for k in range(7, 7 - required_steps, -1): - # Grab 4 bits of each of I and J - offset = k * _S2_LOOKUP_BITS - bits += ((ij[0] >> offset) & lookup_mask) << (_S2_LOOKUP_BITS + 2) - bits += ((ij[1] >> offset) & lookup_mask) << 2 - - # Map bits from iiiijjjjoo to ppppppppoo using lookup table - bits = _S2_LOOKUP_POS[bits] - - # Insert position bits into cell ID - cell_id |= (bits >> 2) << (k * 2 * _S2_LOOKUP_BITS) - - # Remove position bits, leaving just new swap and invert bits for the next round - bits &= _S2_SWAP_MASK | _S2_INVERT_MASK # Mask: 0b11 - - # Left shift and add trailing bit - # The trailing bit addition is disabled, as we are overwriting this below in the truncation - # anyway. This line is kept as an example of the full method for S2 cell ID creation as is done - # in the standard library versions. - cell_id = cell_id << 1 # + 1 - - # Truncate to desired level - # This is done by finding the mask of the trailing 1 bit for the specified level, then zeroing - # out all bits less significant than this, then finally setting the trailing 1 bit. This is - # still necessary to do even after a reduced number of steps `required_steps` above, since each - # step contains multiple levels that may need partial overwrite. Additionally, we need to add - # the trailing 1 bit, which is not yet set above. - least_significant_bit_mask = 1 << (2 * (_S2_MAX_LEVEL - level)) - cell_id = (cell_id & -least_significant_bit_mask) | least_significant_bit_mask - - return cell_id - - -def lat_lon_to_token(lat: float, lon: float, level: int = 30) -> str: - """ - Convert lat/lon to a S2 token. - - Converts the S2 cell ID to hex and strips any trailing zeros. The 0 cell ID token is represented - as 'X' to prevent it being an empty string. - - It is expected that the lat/lon provided are normalised, with latitude in the range -90 to 90. - - See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.cc#L204-L220 - - Args: - lat: The latitude to convert, in degrees. - lon: The longitude to convert, in degrees. - level: The level of the cell ID to generate, from 0 up to 30. - - Returns: - The S2 token string for the lat/lon location. - - Raises: - ValueError: When level is not an integer, is < 0 or is > 30. - - """ - # Generate cell ID and convert to token - return cell_id_to_token(lat_lon_to_cell_id(lat=lat, lon=lon, level=level)) - - -# -# Decode functions -# - -def cell_id_to_lat_lon( # pylint: disable=too-many-locals - cell_id: int -) -> Tuple[float, float]: - """ - Convert S2 cell ID to lat/lon. - - Args: - cell_id: The S2 cell ID integer. - - Returns: - The lat/lon (in degrees) tuple generated from the S2 cell ID. - - Raises: - TypeError: If the cell_id is not int. - InvalidCellID: If the cell_id is invalid. - - """ - # Check input - if not cell_id_is_valid(cell_id): - raise InvalidCellID('Cannot decode invalid S2 cell ID: {}'.format(cell_id)) - - # Populate _S2_LOOKUP_IJ on first run. - # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.cc#L75-L109 - # This table takes 10 bits of curve position and orientation and returns 10 bits of I and J and - # new orientation - if _S2_LOOKUP_IJ is None: # pragma: no cover - _s2_init_lookups() - - # Extract face + IJ from cell ID - # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.cc#L312-L367 - # - # This is done by looking up 8 bits of curve position at a time in the lookup table, along with - # two bits of orientation (swap (1) and invert (2)). This gives back 8 bits of I and J (4 each) - # and two new orientation bits for the curve within the sub-cells in the next step. - # - # The swap bit swaps I and J with each other - # The invert bit inverts the bits of I and J, which means axes are negated - # - # In the first loop (most significant bits), the 3 bits occupied by the face need to be masked - # out, since these are not set in the IJ to cell ID during encoding. - # - # The I and J returned here are of one of the two leaf (level 30) cells that are located - # diagonally closest to the cell center. This happens because repeated ..00.. will select the - # 'lower left' (for nominally oriented Hilbert curve segments) of the sub-cells. The ..10.. - # arising from the trailing bit, prior to the repeated ..00.. bits, ensures we first pick the - # 'upper right' of the cell, then iterate in to lower left until we hit the leaf cell. This - # means we pick the leaf cell to the north east of the parent cell center (again for nominal - # orientation). - # However, in the case of the swapped and inverted curve segment (4th sub-curve segment), the - # ..10.. will select the 'lower left' and then iterate to the 'upper right' with each ..00.. - # following. In that case, we will be offset left and down by one leaf cell in each of I and J, - # which needs to be fixed to have a consistent mapping. This is detectable by seeing that the - # final bit of I or J is 1 (i.e we have picked an odd row/column, which will happen concurrently - # in both I and J, so we only need to check one), except in case of level 29 where the logic is - # inverted and the correction needs to be applied when we pick an even row/column (i.e I/J ends - # in 0), since there are no trailing ..00.. available after the ``..10..`` when we are at level - # 29+. - # - # This behaviour can be captured in the expression: - # apply_correction = not leaf and (i ^ (is level 29)) & 1 - # apply_correction = not leaf and (i ^ (cell_id >> 2)) & 1 - # - # We check for level 29 by looking for the trailing 1 in the third LSB, when we already know - # that we are not a leaf cell (which could give false positive) by the initial check in the - # expression. - # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.h#L503-L529 - # - face = cell_id >> _S2_POS_BITS - bits = face & _S2_SWAP_MASK # ppppppppoo. Initially set by by face - lookup_mask = (1 << _S2_LOOKUP_BITS) - 1 # Mask of 4 one bits: 0b1111 - i = 0 - j = 0 - for k in range(7, -1, -1): - # Pull out 8 bits of cell ID, except in first loop where we pull out only 4 - n_bits = (_S2_MAX_LEVEL - 7 * _S2_LOOKUP_BITS) if k == 7 else _S2_LOOKUP_BITS - extract_mask = (1 << (2 * n_bits)) - 1 # 8 (or 4) one bits - bits += ( - (cell_id >> (k * 2 * _S2_LOOKUP_BITS + 1)) & extract_mask - ) << 2 - - # Map bits from ppppppppoo to iiiijjjjoo using lookup table - bits = _S2_LOOKUP_IJ[bits] - - # Extract I and J bits - offset = k * _S2_LOOKUP_BITS - i += (bits >> (_S2_LOOKUP_BITS + 2)) << offset # Don't need lookup mask here - j += ((bits >> 2) & lookup_mask) << offset - - # Remove I and J bits, leaving just new swap and invert bits for the next round - bits &= _S2_SWAP_MASK | _S2_INVERT_MASK # Mask: 0b11 - - # Resolve the center of the cell. For leaf cells, we add half the leaf cell size. For non-leaf - # cells, we currently have one of either two cells diagonally around the cell center and want - # to pick the leaf-cell edges that represent the parent cell center, as described above. The - # center_correction_delta is 2x the offset, as we left shift I and J first. - # This gives us the values Si and Ti, which are discrete representation of S and T in range 0 to - # _S2_MAX_SI_TI. The extra power of 2 over IJ allows for identifying both the center and edge of - # cells, whilst IJ is just the leaf cells. - # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2coords.h#L57-L65 - is_leaf = bool(cell_id & 1) # Cell is leaf cell when trailing one bit is in LSB - apply_correction = not is_leaf and ((i ^ (cell_id >> 2)) & 1) - correction_delta = 1 if is_leaf else (2 if apply_correction else 0) - si = (i << 1) + correction_delta # pylint: disable=invalid-name - ti = (j << 1) + correction_delta # pylint: disable=invalid-name - - # Convert integer si/ti to double ST - # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2coords.h#L338-L341 - st = (_s2_si_ti_to_st(si), _s2_si_ti_to_st(ti)) # pylint: disable=invalid-name - - # Project cell-space ST to cube-space UV - # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2coords.h#L312-L315 - uv = (_s2_st_to_uv(st[0]), _s2_st_to_uv(st[1])) # pylint: disable=invalid-name - - # Convert face + UV to S2Point XYZ - # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2coords.h#L348-L357 - s2_point = _s2_face_uv_to_xyz(face, uv) - - # Normalise XYZ S2Point vector - # This section is part of the reference implementation but is not necessary when mapping - # straight into lat/lon, since the normalised and unnormalised triangles used to calculate the - # angles are geometrically similar. If anything, the normalisation process loses precision when - # tested against the reference implementation, albeit not at a level that is important either - # way. The code below is left for demonstration of the normalisation process. - # norm = math.sqrt(s2_point[0] ** 2 + s2_point[1] ** 2 + s2_point[2] ** 2) - # s2_point = (s2_point[0] / norm, s2_point[1] / norm, s2_point[2] / norm) - - # Map into lat/lon - # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2latlng.h#L196-L205 - lat_rad = math.atan2(s2_point[2], math.sqrt(s2_point[0] ** 2 + s2_point[1] ** 2)) - lon_rad = math.atan2(s2_point[1], s2_point[0]) - - return (math.degrees(lat_rad), math.degrees(lon_rad)) - - -def token_to_lat_lon(token: str) -> Tuple[float, float]: - """ - Convert S2 token to lat/lon. - - See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.cc#L222-L239 - - Args: - token: The S2 token string. Can be upper or lower case hex string. - - Returns: - The lat/lon (in degrees) tuple generated from the S2 token. - - Raises: - TypeError: If the token is not str. - InvalidToken: If the token length is over 16. - InvalidToken: If the token is invalid. - InvalidCellID: If the contained cell_id is invalid. - - """ - # Check input - if not token_is_valid(token): - raise InvalidToken('Cannot decode invalid S2 token: {}'.format(token)) - - # Convert to cell ID and decode to lat/lon - return cell_id_to_lat_lon(token_to_cell_id(token)) - - -# -# Token canonicalisation -# - -def token_to_canonical_token(token: str) -> str: - """ - Convert S2 token to a canonicalised S2 token. - - This produces a token that matches the form generated by the reference C++ implementation: - - - Lower case (except 'X' below) - - No whitespace - - Trailing '0' characters stripped - - Zero cell ID represented as 'X', not 'x' or '' - - Args: - token: The S2 token string to canonicalise. - - Returns: - The canonicalised S2 token. - - """ - # Convert token to lower case. - # Note that 'X' below will be returned upper case - token = token.lower() - - # Strip any surrounding whitespace - token = token.strip() - - # Strip any trailing zeros - token = token.rstrip('0') - - # If empty string or 'x', return 'X' token - if token in ('', 'x'): - token = 'X' - - return token - - -# -# Validation -# - -def cell_id_is_valid(cell_id: int) -> bool: - """ - Check that a S2 cell ID is valid. - - Looks for valid face bits and a trailing 1 bit in one of the correct locations. - - Args: - cell_id: The S2 cell integer to validate. - - Returns: - True if the cell ID is valid, False otherwise. - - Raises: - TypeError: If the cell_id is not int. - - """ - # Check input - if not isinstance(cell_id, int): - raise TypeError('Cannot decode S2 cell ID from type: {}'.format(type(cell_id))) - - # Check for zero ID - # This avoids overflow warnings below when 1 gets added to max uint64 - if cell_id == 0: - return False - - # Check face bits - if (cell_id >> _S2_POS_BITS) > 5: - return False - - # Check trailing 1 bit is in one of the even bit positions allowed for the 30 levels, using the - # mask: 0b0001010101010101010101010101010101010101010101010101010101010101 = 0x1555555555555555 - lowest_set_bit = cell_id & (~cell_id + 1) # pylint: disable=invalid-unary-operand-type - if not lowest_set_bit & 0x1555555555555555: - return False - - return True # Checks have passed, cell ID must be valid - - -def token_is_valid(token: str) -> bool: - """ - Check that a S2 token is valid. - - Looks for valid characters, then checks that the contained S2 cell ID is also valid. Note that - the '', 'x' and 'X' tokens are considered invalid, since the cell IDs they represent are - invalid. - - Args: - token: The S2 token string to validate. - - Returns: - True if the token is valid, False otherwise. - - Raises: - TypeError: If the token is not str. - - """ - # Check input - if not isinstance(token, str): - raise TypeError('Cannot check S2 token with type: {}'.format(type(token))) - - # First check string with regex - if not re.match(r'^[0-9a-fA-f]{1,16}$', token): - return False - - # Check the contained cell ID is valid - return cell_id_is_valid(token_to_cell_id(token)) - - -# -# Level extraction functions -# - -def cell_id_to_level(cell_id: int) -> int: - """ - Get the level for a S2 cell ID. - - See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.h#L543-L551 - - Args: - cell_id: The S2 cell ID integer. - - Returns: - The level of the S2 cell ID. - - Raises: - TypeError: If the cell_id is not int. - InvalidCellID: If the cell_id is invalid. - - """ - # Check input - if not cell_id_is_valid(cell_id): - raise InvalidCellID('Cannot decode invalid S2 cell ID: {}'.format(cell_id)) - - # Find the position of the lowest set one bit, which will be the trailing one bit. The level is - # given by the max level (30) minus the floored division by two of the position of the lowest - # set bit. - # - # The position of the lowest set bit is found using 'count trailing zeros', which would be - # equivalent to the C++20 function std::countr_zero() or the ctz instruction. - lsb_pos = 0 - while cell_id != 0: # pragma: no branch - if cell_id & 1: - break - lsb_pos += 1 - cell_id >>= 1 - - return int(_S2_MAX_LEVEL - (lsb_pos >> 1)) - - -def token_to_level(token: str) -> int: - """ - Get the level for a S2 token. - - See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.h#L543-L551 - - Args: - token: The S2 token string. Can be upper or lower case hex string. - - Returns: - The level of the S2 token. - - Raises: - TypeError: If the token is not str. - InvalidToken: If the token length is over 16. - InvalidToken: If the token is invalid. - InvalidCellID: If the contained cell_id is invalid. - - """ - # Check input - if not token_is_valid(token): - raise InvalidToken('Cannot decode invalid S2 token: {}'.format(token)) - - # Convert to cell ID and get the level for that - return cell_id_to_level(token_to_cell_id(token)) - - -# -# Parent functions -# - -def cell_id_to_parent_cell_id( - cell_id: int, level: Optional[int] = None -) -> int: - """ - Get the parent cell ID of a S2 cell ID. - - Args: - cell_id: The S2 cell ID integer. - level: The parent level to get the cell ID for. Must be less than or equal to the current - level of the provided cell ID. If unspecified, or None, the direct parent cell ID will - be returned. - - Returns: - The parent cell ID at the specified level. - - Raises: - TypeError: If the cell_id is not int. - InvalidCellID: If the cell_id is invalid. - ValueError: If cell ID is already level 0 and level is None. - ValueError: When level is not an integer, is < 0 or is > 30. - ValueError: If level is greater than the provided cell ID level. - - """ - # Check input - if not cell_id_is_valid(cell_id): - raise InvalidCellID('Cannot decode invalid S2 cell ID: {}'.format(cell_id)) - - # Get current level of the cell ID and check it is suitable with the requested level - current_level = cell_id_to_level(cell_id) - if level is None and current_level == 0: - raise ValueError('Cannot get parent cell ID of a level 0 cell ID') - if level is None: - level = current_level - 1 - - if not isinstance(level, int) or level < 0 or level > _S2_MAX_LEVEL: - raise ValueError('S2 level must be integer >= 0 and <= 30') - - if level > current_level: - raise ValueError('Cannot get level {} parent cell ID of cell ID with level {}'.format( - level, current_level - )) - if level == current_level: - # Requested parent level is current level, return cell ID itself - return cell_id - - # Truncate to desired level - # This is done by finding the mask of the trailing 1 bit for the specified level, then zeroing - # out all bits less significant than this, then finally setting the trailing 1 bit. This is - # still necessary to do even after a reduced number of steps `required_steps` above, since each - # step contains multiple levels that may need partial overwrite. Additionally, we need to add - # the trailing 1 bit, which is not yet set above. - least_significant_bit_mask = 1 << (2 * (_S2_MAX_LEVEL - level)) - cell_id = (cell_id & -least_significant_bit_mask) | least_significant_bit_mask - - return cell_id - - -def token_to_parent_token(token: str, level: Optional[int] = None) -> str: - """ - Get the parent token of a S2 token. - - Args: - token: The S2 token string. Can be upper or lower case hex string. - level: The parent level to get the token for. Must be less than or equal to the current - level of the provided toke. If unspecified, or None, the direct parent token will be - returned. - - Returns: - The parent token at the specified level. - - Raises: - TypeError: If the token is not str. - InvalidToken: If the token length is over 16. - InvalidToken: If the token is invalid. - InvalidCellID: If the contained cell_id is invalid. - ValueError: If token is already level 0 and level is None. - ValueError: When level is not an integer, is < 0 or is > 30. - ValueError: If level is greater than the provided token level. - - """ - # Check input - if not token_is_valid(token): - raise InvalidToken('Cannot decode invalid S2 token: {}'.format(token)) - - # Convert to cell ID and get parent and convert back to token - return cell_id_to_token(cell_id_to_parent_cell_id(token_to_cell_id(token), level)) +__version__ = '1.5.0' diff --git a/s2cell/s2cell.py b/s2cell/s2cell.py new file mode 100644 index 0000000..575b301 --- /dev/null +++ b/s2cell/s2cell.py @@ -0,0 +1,947 @@ +# Copyright 2020 Adam Liddell +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Minimal Python S2 Geometry cell ID, token and lat/lon conversion library.""" + +import math +import re +from typing import Optional, Tuple + + +# +# s2cell exceptions +# + +class InvalidCellID(Exception): + """Exception type for invalid cell IDs.""" + + +class InvalidToken(Exception): + """Exception type for invalid tokens.""" + + +# +# S2 base constants needed for cell mapping +# + +# The maximum level supported within an S2 cell ID. Each level is represented by two bits in the +# final cell ID +_S2_MAX_LEVEL = 30 + +# The maximum value within the I and J bits of an S2 cell ID +_S2_MAX_SIZE = 1 << _S2_MAX_LEVEL + +# The number of bits in a S2 cell ID used for specifying the base face +_S2_FACE_BITS = 3 + +# The number of bits in a S2 cell ID used for specifying the position along the Hilbert curve +_S2_POS_BITS = 2 * _S2_MAX_LEVEL + 1 + +# The maximum value of the Si/Ti integers used when mapping from IJ to ST. This is twice the max +# value of I and J, since Si/Ti allow referencing both the center and edge of a leaf cell +_S2_MAX_SI_TI = 1 << (_S2_MAX_LEVEL + 1) + +# Mask that specifies the swap orientation bit for the Hilbert curve +_S2_SWAP_MASK = 1 + +# Mask that specifies the invert orientation bit for the Hilbert curve +_S2_INVERT_MASK = 2 + +# The number of bits per I and J in the lookup tables +_S2_LOOKUP_BITS = 4 + +# Lookup table for mapping 10 bits of IJ + orientation to 10 bits of Hilbert curve position + +# orientation. Populated later by _s2_init_lookups +_S2_LOOKUP_POS = None + +# Lookup table for mapping 10 bits of Hilbert curve position + orientation to 10 bits of IJ + +# orientation. Populated later by _s2_init_lookups +_S2_LOOKUP_IJ = None + +# Lookup table of two bits of IJ from two bits of curve position, based also on the current curve +# orientation from the swap and invert bits +_S2_POS_TO_IJ = [ + [0, 1, 3, 2], # 0: Normal order, no swap or invert + [0, 2, 3, 1], # 1: Swap bit set, swap I and J bits + [3, 2, 0, 1], # 2: Invert bit set, invert bits + [3, 1, 0, 2], # 3: Swap and invert bits set +] + +# Lookup for the orientation update mask of one of the four sub-cells within a higher level cell. +# This mask is XOR'ed with the current orientation to get the sub-cell orientation. +_S2_POS_TO_ORIENTATION_MASK = [_S2_SWAP_MASK, 0, 0, _S2_SWAP_MASK | _S2_INVERT_MASK] + + +# +# S2 helper functions +# + +def _s2_uv_to_st(component: float) -> float: + """ + Convert S2 UV to ST. + + This is done using the quadratic projection that is used by default for S2. The C++ and Java S2 + libraries use a different definition of the ST cell-space, but the end result in IJ is the same. + The below uses the C++ ST definition. + + See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2coords.h#L317-L320 + + """ + if component >= 0.0: + return 0.5 * math.sqrt(1.0 + 3.0 * component) + return 1.0 - 0.5 * math.sqrt(1.0 - 3.0 * component) + + +def _s2_st_to_uv(component: float) -> float: + """ + Convert S2 ST to UV. + + This is done using the quadratic projection that is used by default for S2. The C++ and Java S2 + libraries use a different definition of the ST cell-space, but the end result in IJ is the same. + The below uses the C++ ST definition. + + See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2coords.h#L312-L315 + + """ + if component >= 0.5: + return (1.0 / 3.0) * (4.0 * component ** 2 - 1.0) + return (1.0 / 3.0) * (1.0 - 4.0 * (1.0 - component) ** 2) + + +def _s2_st_to_ij(component: float) -> int: + """ + Convert S2 ST to IJ. + + The mapping here differs between C++ and Java versions, but the combination of + _st_to_ij(_uv_to_st(val)) is the same for both. The below uses the C++ ST definition. + + See s2geometry/blob/2c02e21040e0b82aa5719e96033d02b8ce7c0eff/src/s2/s2coords.h#L333-L336 + + """ + # The reference implementation does round(_S2_MAX_SIZE * component - 0.5), which is equivalent + # to math.floor(_S2_MAX_SIZE * component) + return int(max(0, min(_S2_MAX_SIZE - 1, math.floor(_S2_MAX_SIZE * component)))) + + +def _s2_si_ti_to_st(component: int) -> float: + """ + Convert S2 Si/Ti to ST. + + This converts an integer in range 0 to _S2_MAX_SI_TI into a float in range 0.0 to 1.0. + + See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2coords.h#L338-L341 + + """ + return (1.0 / _S2_MAX_SI_TI) * component + + +def _s2_face_uv_to_xyz( # pylint: disable=invalid-name + face: int, uv: Tuple[float, float] +) -> Tuple[float, float, float]: + """ + Convert face + UV to S2Point XYZ. + + See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2coords.h#L348-L357 + + Args: + face: The S2 face for the input point. + uv: The S2 face UV coordinates. + + Returns: + The unnormalised S2Point XYZ. + + Raises: + ValueError: If the face is not valid in range 0-5. + + """ + # Face -> XYZ components -> indices with negation: + # 0 -> ( 1, u, v) -> ( /, 0, 1) + # 1 -> (-u, 1, v) -> (-0, /, 1) + # 2 -> (-u, -v, 1) -> (-0, -1, /) + # 3 -> (-1, -v, -u) -> (-/, -1, -0) <- -1 here means -1 times the value in index 1, + # 4 -> ( v, -1, -u) -> ( 1, -/, -0) not index -1 + # 5 -> ( v, u, -1) -> ( 1, 0, -/) + if face == 0: + s2_point = (1, uv[0], uv[1]) + elif face == 1: + s2_point = (-uv[0], 1, uv[1]) + elif face == 2: + s2_point = (-uv[0], -uv[1], 1) + elif face == 3: + s2_point = (-1, -uv[1], -uv[0]) + elif face == 4: + s2_point = (uv[1], -1, -uv[0]) + elif face == 5: + s2_point = (uv[1], uv[0], -1) + else: + raise ValueError('Cannot convert UV to XYZ with invalid face: {}'.format(face)) + + return s2_point + + +def _s2_init_lookups() -> None: + """ + Initialise the S2 lookups in global vars _S2_LOOKUP_POS and _S2_LOOKUP_IJ. + + This generates 4 variations of a 4 level deep Hilbert curve, one for each swap/invert bit + combination. This allows mapping between 8 bits (+2 orientation) of Hilbert curve position and 8 + bits (+2 orientation) of I and J, and vice versa. The new orientation bits read from the mapping + tell us the base orientation of the curve segments within the next deeper level of sub-cells. + + This implementation differs in structure from the reference implementation, since it is + iterative rather than recursive. The end result is the same lookup table. + + See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.cc#L75-L109 + + """ + global _S2_LOOKUP_POS, _S2_LOOKUP_IJ # pylint: disable=global-statement + if _S2_LOOKUP_POS is None or _S2_LOOKUP_IJ is None: # pragma: no branch + # Initialise empty lookup tables + lookup_length = 1 << int(2 * _S2_LOOKUP_BITS + 2) # = 1024 + _S2_LOOKUP_POS = [0] * lookup_length + _S2_LOOKUP_IJ = [0] * lookup_length + + # Generate lookups for each of the base orientations given by the swap and invert bits + for base_orientation in [ + 0, _S2_SWAP_MASK, _S2_INVERT_MASK, _S2_SWAP_MASK | _S2_INVERT_MASK # 0-3 effectively + ]: + # Walk the 256 possible positions within a level 4 curve. This implementation is not + # the fastest since it does not reuse the common ancestor of neighbouring positions, but + # is simpler to read + for pos in range(4 ** 4): # 4 levels of sub-divisions + ij = 0 # Has pattern iiiijjjj, not ijijijij # pylint: disable=invalid-name + orientation = base_orientation + + # Walk the pairs of bits of pos, from most significant to least, getting IJ and + # orientation as we go + for bit_pair_offset in range(4): + # Bit pair is effectively the sub-cell index + bit_pair = (pos >> ((3 - bit_pair_offset) * 2)) & 0b11 + + # Get the I and J for the sub-cell index. These need to be spread into iiiijjjj + # by inserting as bit positions 4 and 0 + ij_bits = _S2_POS_TO_IJ[orientation][bit_pair] + ij = ( # pylint: disable=invalid-name + (ij << 1) # Free up position 4 and 0 from old IJ + | ((ij_bits & 2) << 3) # I bit in position 4 + | (ij_bits & 1) # J bit in position 0 + ) + + # Update the orientation with the new sub-cell orientation + orientation = orientation ^ _S2_POS_TO_ORIENTATION_MASK[bit_pair] + + # Shift IJ and position to allow orientation bits in LSBs of lookup + ij <<= 2 # pylint: disable=invalid-name + pos <<= 2 + + # Write lookups + _S2_LOOKUP_POS[ij | base_orientation] = pos | orientation + _S2_LOOKUP_IJ[pos | base_orientation] = ij | orientation + + +# +# Cell ID <-> Token translation functions +# + +def cell_id_to_token(cell_id: int) -> str: + """ + Convert S2 cell ID to a S2 token. + + Converts the S2 cell ID to hex and strips any trailing zeros. The 0 cell ID token is represented + as 'X' to prevent it being an empty string. + + See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.cc#L204-L220 + + Args: + cell_id: The S2 cell ID integer. + + Returns: + The S2 token string for the S2 cell ID. + + Raises: + TypeError: If the cell_id is not int. + + """ + # Check type + if not isinstance(cell_id, int): + raise TypeError('Cannot convert S2 cell ID from type: {}'.format(type(cell_id))) + + # The zero token is encoded as 'X' rather than as a zero-length string + if cell_id == 0: + return 'X' + + # Convert cell ID to 16 character hex string and strip any implicit trailing zeros + return '{:016x}'.format(cell_id).rstrip('0') + + +def token_to_cell_id(token: str) -> int: + """ + Convert S2 token to S2 cell ID. + + Restores the stripped 0 characters from the token and converts the hex string to integer. + + See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.cc#L222-L239 + + Args: + token: The S2 token string. Can be upper or lower case hex string. + + Returns: + The S2 cell ID for the S2 token. + + Raises: + TypeError: If the token is not str. + InvalidToken: If the token length is over 16. + + """ + # Check input + if not isinstance(token, str): + raise TypeError('Cannot convert S2 token from type: {}'.format(type(token))) + + if len(token) > 16: + raise InvalidToken('Cannot convert S2 token with length > 16 characters') + + # Check for the zero cell ID represented by the character 'x' or 'X' rather than as the empty + # string + if token in ('x', 'X'): + return 0 + + # Add stripped implicit zeros to create the full 16 character hex string + token = token + ('0' * (16 - len(token))) + + # Convert to cell ID by converting hex to int + return int(token, 16) + + +# +# Encode functions +# + +def lat_lon_to_cell_id( # pylint: disable=too-many-locals + lat: float, lon: float, level: int = 30 +) -> int: + """ + Convert lat/lon to a S2 cell ID. + + It is expected that the lat/lon provided are normalised, with latitude in the range -90 to 90. + + Args: + lat: The latitude to convert, in degrees. + lon: The longitude to convert, in degrees. + level: The level of the cell ID to generate, from 0 up to 30. + + Returns: + The S2 cell ID for the lat/lon location. + + Raises: + ValueError: When level is not an integer, is < 0 or is > 30. + + """ + if not isinstance(level, int) or level < 0 or level > _S2_MAX_LEVEL: + raise ValueError('S2 level must be integer >= 0 and <= 30') + + # Populate _S2_LOOKUP_POS on first run. + # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.cc#L75-L109 + # + # This table takes 10 bits of I and J and orientation and returns 10 bits of curve position and + # new orientation + if _S2_LOOKUP_POS is None: # pragma: no cover + _s2_init_lookups() + + # Reuse constant expressions + lat_rad = math.radians(lat) + lon_rad = math.radians(lon) + sin_lat_rad = math.sin(lat_rad) + cos_lat_rad = math.cos(lat_rad) + sin_lon_rad = math.sin(lon_rad) + cos_lon_rad = math.cos(lon_rad) + + # Convert to S2Point + # This is effectively the unit non-geodetic ECEF vector + s2_point = ( + cos_lat_rad * cos_lon_rad, # X + cos_lat_rad * sin_lon_rad, # Y + sin_lat_rad # Z + ) + + # Get cube face + # See s2geometry/blob/2c02e21040e0b82aa5719e96033d02b8ce7c0eff/src/s2/s2coords.h#L380-L384 + # + # The face is determined by the largest XYZ component of the S2Point vector. When the component + # is negative, the second set of three faces is used. + # Largest component -> face: + # +x -> 0 + # +y -> 1 + # +z -> 2 + # -x -> 3 + # -y -> 4 + # -z -> 5 + face = max(enumerate(s2_point), key=lambda p: abs(p[1]))[0] # Largest absolute component + if s2_point[face] < 0.0: + face += 3 + + # Convert face + XYZ to cube-space face + UV + # See s2geometry/blob/2c02e21040e0b82aa5719e96033d02b8ce7c0eff/src/s2/s2coords.h#L366-L372 + # + # The faces are oriented to ensure continuity of curve. + # Face -> UV components -> indices with negation (without divisor, which is always the remaining + # component (index: face % 3)): + # 0 -> ( y, z) -> ( 1, 2) + # 1 -> (-x, z) -> (-0, 2) + # 2 -> (-x, -y) -> (-0, -1) <- -1 here means -1 times the value in index 1, not index -1 + # 3 -> ( z, y) -> ( 2, 1) + # 4 -> ( z, -x) -> ( 2, -0) + # 5 -> (-y, -x) -> (-1, -0) + # + # For a compiled language, a switch statement on face is preferable as it will be more easily + # optimised as a jump table etc; but in Python the indexing method is more concise. + # + # The index selection can be reduced to some bit magic: + # U: 1 - ((face + 1) >> 1) + # V: 2 - (face >> 1) + # + # The negation of the the two components is then selected: + # U: (face in [1, 2, 5]) ? -1: 1 + # V: (face in [2, 4, 5])) ? -1: 1 + uv = ( # pylint: disable=invalid-name + s2_point[1 - ((face + 1) >> 1)] / s2_point[face % 3], # U + s2_point[2 - (face >> 1)] / s2_point[face % 3] # V + ) + if face in (1, 2, 5): + uv = (-uv[0], uv[1]) # Negate U # pylint: disable=invalid-name + if face in (2, 4, 5): + uv = (uv[0], -uv[1]) # Negate V # pylint: disable=invalid-name + + # Project cube-space UV to cell-space ST + # See s2geometry/blob/2c02e21040e0b82aa5719e96033d02b8ce7c0eff/src/s2/s2coords.h#L317-L320 + st = (_s2_uv_to_st(uv[0]), _s2_uv_to_st(uv[1])) # pylint: disable=invalid-name + + # Convert ST to IJ integers + # See s2geometry/blob/2c02e21040e0b82aa5719e96033d02b8ce7c0eff/src/s2/s2coords.h#L333-L336 + ij = (_s2_st_to_ij(st[0]), _s2_st_to_ij(st[1])) # pylint: disable=invalid-name + + # Convert face + IJ to cell ID + # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.cc#L256-L298 + # + # This is done by looking up 8 bits of I and J (4 each) at a time in the lookup table, along + # with two bits of orientation (swap (1) and invert (2)). This gives back 8 bits of position + # along the curve and two new orientation bits for the curve within the sub-cells in the next + # step. + # + # The swap bit swaps I and J with each other + # The invert bit inverts the bits of I and J, which means axes are negated + # + # Compared to the standard versions, we check the required number of steps we need to do for the + # requested level and don't perform steps that will be completely overwritten in the truncation + # below, rather than always doing every step. Each step does 4 bits each of I and J, which is 4 + # levels, so the required number of steps is ceil((level + 2) / 4), when level is > 0. The + # additional 2 levels added are required to account for the top 3 bits (4 before right shift) + # that are occupied by the face bits. + bits = face & _S2_SWAP_MASK # iiiijjjjoo. Initially set by by face + cell_id = face << (_S2_POS_BITS - 1) # Insert face at second most signficant bits + lookup_mask = (1 << int(_S2_LOOKUP_BITS)) - 1 # Mask of 4 one bits: 0b1111 + required_steps = math.ceil((level + 2) / 4) if level > 0 else 0 + for k in range(7, 7 - required_steps, -1): + # Grab 4 bits of each of I and J + offset = k * _S2_LOOKUP_BITS + bits += ((ij[0] >> offset) & lookup_mask) << (_S2_LOOKUP_BITS + 2) + bits += ((ij[1] >> offset) & lookup_mask) << 2 + + # Map bits from iiiijjjjoo to ppppppppoo using lookup table + bits = _S2_LOOKUP_POS[bits] + + # Insert position bits into cell ID + cell_id |= (bits >> 2) << (k * 2 * _S2_LOOKUP_BITS) + + # Remove position bits, leaving just new swap and invert bits for the next round + bits &= _S2_SWAP_MASK | _S2_INVERT_MASK # Mask: 0b11 + + # Left shift and add trailing bit + # The trailing bit addition is disabled, as we are overwriting this below in the truncation + # anyway. This line is kept as an example of the full method for S2 cell ID creation as is done + # in the standard library versions. + cell_id = cell_id << 1 # + 1 + + # Truncate to desired level + # This is done by finding the mask of the trailing 1 bit for the specified level, then zeroing + # out all bits less significant than this, then finally setting the trailing 1 bit. This is + # still necessary to do even after a reduced number of steps `required_steps` above, since each + # step contains multiple levels that may need partial overwrite. Additionally, we need to add + # the trailing 1 bit, which is not yet set above. + least_significant_bit_mask = 1 << (2 * (_S2_MAX_LEVEL - level)) + cell_id = (cell_id & -least_significant_bit_mask) | least_significant_bit_mask + + return cell_id + + +def lat_lon_to_token(lat: float, lon: float, level: int = 30) -> str: + """ + Convert lat/lon to a S2 token. + + Converts the S2 cell ID to hex and strips any trailing zeros. The 0 cell ID token is represented + as 'X' to prevent it being an empty string. + + It is expected that the lat/lon provided are normalised, with latitude in the range -90 to 90. + + See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.cc#L204-L220 + + Args: + lat: The latitude to convert, in degrees. + lon: The longitude to convert, in degrees. + level: The level of the cell ID to generate, from 0 up to 30. + + Returns: + The S2 token string for the lat/lon location. + + Raises: + ValueError: When level is not an integer, is < 0 or is > 30. + + """ + # Generate cell ID and convert to token + return cell_id_to_token(lat_lon_to_cell_id(lat=lat, lon=lon, level=level)) + + +# +# Decode functions +# + +def cell_id_to_lat_lon( # pylint: disable=too-many-locals + cell_id: int +) -> Tuple[float, float]: + """ + Convert S2 cell ID to lat/lon. + + Args: + cell_id: The S2 cell ID integer. + + Returns: + The lat/lon (in degrees) tuple generated from the S2 cell ID. + + Raises: + TypeError: If the cell_id is not int. + InvalidCellID: If the cell_id is invalid. + + """ + # Check input + if not cell_id_is_valid(cell_id): + raise InvalidCellID('Cannot decode invalid S2 cell ID: {}'.format(cell_id)) + + # Populate _S2_LOOKUP_IJ on first run. + # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.cc#L75-L109 + # This table takes 10 bits of curve position and orientation and returns 10 bits of I and J and + # new orientation + if _S2_LOOKUP_IJ is None: # pragma: no cover + _s2_init_lookups() + + # Extract face + IJ from cell ID + # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.cc#L312-L367 + # + # This is done by looking up 8 bits of curve position at a time in the lookup table, along with + # two bits of orientation (swap (1) and invert (2)). This gives back 8 bits of I and J (4 each) + # and two new orientation bits for the curve within the sub-cells in the next step. + # + # The swap bit swaps I and J with each other + # The invert bit inverts the bits of I and J, which means axes are negated + # + # In the first loop (most significant bits), the 3 bits occupied by the face need to be masked + # out, since these are not set in the IJ to cell ID during encoding. + # + # The I and J returned here are of one of the two leaf (level 30) cells that are located + # diagonally closest to the cell center. This happens because repeated ..00.. will select the + # 'lower left' (for nominally oriented Hilbert curve segments) of the sub-cells. The ..10.. + # arising from the trailing bit, prior to the repeated ..00.. bits, ensures we first pick the + # 'upper right' of the cell, then iterate in to lower left until we hit the leaf cell. This + # means we pick the leaf cell to the north east of the parent cell center (again for nominal + # orientation). + # However, in the case of the swapped and inverted curve segment (4th sub-curve segment), the + # ..10.. will select the 'lower left' and then iterate to the 'upper right' with each ..00.. + # following. In that case, we will be offset left and down by one leaf cell in each of I and J, + # which needs to be fixed to have a consistent mapping. This is detectable by seeing that the + # final bit of I or J is 1 (i.e we have picked an odd row/column, which will happen concurrently + # in both I and J, so we only need to check one), except in case of level 29 where the logic is + # inverted and the correction needs to be applied when we pick an even row/column (i.e I/J ends + # in 0), since there are no trailing ..00.. available after the ``..10..`` when we are at level + # 29+. + # + # This behaviour can be captured in the expression: + # apply_correction = not leaf and (i ^ (is level 29)) & 1 + # apply_correction = not leaf and (i ^ (cell_id >> 2)) & 1 + # + # We check for level 29 by looking for the trailing 1 in the third LSB, when we already know + # that we are not a leaf cell (which could give false positive) by the initial check in the + # expression. + # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.h#L503-L529 + # + face = cell_id >> _S2_POS_BITS + bits = face & _S2_SWAP_MASK # ppppppppoo. Initially set by by face + lookup_mask = (1 << _S2_LOOKUP_BITS) - 1 # Mask of 4 one bits: 0b1111 + i = 0 + j = 0 + for k in range(7, -1, -1): + # Pull out 8 bits of cell ID, except in first loop where we pull out only 4 + n_bits = (_S2_MAX_LEVEL - 7 * _S2_LOOKUP_BITS) if k == 7 else _S2_LOOKUP_BITS + extract_mask = (1 << (2 * n_bits)) - 1 # 8 (or 4) one bits + bits += ( + (cell_id >> (k * 2 * _S2_LOOKUP_BITS + 1)) & extract_mask + ) << 2 + + # Map bits from ppppppppoo to iiiijjjjoo using lookup table + bits = _S2_LOOKUP_IJ[bits] + + # Extract I and J bits + offset = k * _S2_LOOKUP_BITS + i += (bits >> (_S2_LOOKUP_BITS + 2)) << offset # Don't need lookup mask here + j += ((bits >> 2) & lookup_mask) << offset + + # Remove I and J bits, leaving just new swap and invert bits for the next round + bits &= _S2_SWAP_MASK | _S2_INVERT_MASK # Mask: 0b11 + + # Resolve the center of the cell. For leaf cells, we add half the leaf cell size. For non-leaf + # cells, we currently have one of either two cells diagonally around the cell center and want + # to pick the leaf-cell edges that represent the parent cell center, as described above. The + # center_correction_delta is 2x the offset, as we left shift I and J first. + # This gives us the values Si and Ti, which are discrete representation of S and T in range 0 to + # _S2_MAX_SI_TI. The extra power of 2 over IJ allows for identifying both the center and edge of + # cells, whilst IJ is just the leaf cells. + # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2coords.h#L57-L65 + is_leaf = bool(cell_id & 1) # Cell is leaf cell when trailing one bit is in LSB + apply_correction = not is_leaf and ((i ^ (cell_id >> 2)) & 1) + correction_delta = 1 if is_leaf else (2 if apply_correction else 0) + si = (i << 1) + correction_delta # pylint: disable=invalid-name + ti = (j << 1) + correction_delta # pylint: disable=invalid-name + + # Convert integer si/ti to double ST + # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2coords.h#L338-L341 + st = (_s2_si_ti_to_st(si), _s2_si_ti_to_st(ti)) # pylint: disable=invalid-name + + # Project cell-space ST to cube-space UV + # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2coords.h#L312-L315 + uv = (_s2_st_to_uv(st[0]), _s2_st_to_uv(st[1])) # pylint: disable=invalid-name + + # Convert face + UV to S2Point XYZ + # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2coords.h#L348-L357 + s2_point = _s2_face_uv_to_xyz(face, uv) + + # Normalise XYZ S2Point vector + # This section is part of the reference implementation but is not necessary when mapping + # straight into lat/lon, since the normalised and unnormalised triangles used to calculate the + # angles are geometrically similar. If anything, the normalisation process loses precision when + # tested against the reference implementation, albeit not at a level that is important either + # way. The code below is left for demonstration of the normalisation process. + # norm = math.sqrt(s2_point[0] ** 2 + s2_point[1] ** 2 + s2_point[2] ** 2) + # s2_point = (s2_point[0] / norm, s2_point[1] / norm, s2_point[2] / norm) + + # Map into lat/lon + # See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2latlng.h#L196-L205 + lat_rad = math.atan2(s2_point[2], math.sqrt(s2_point[0] ** 2 + s2_point[1] ** 2)) + lon_rad = math.atan2(s2_point[1], s2_point[0]) + + return (math.degrees(lat_rad), math.degrees(lon_rad)) + + +def token_to_lat_lon(token: str) -> Tuple[float, float]: + """ + Convert S2 token to lat/lon. + + See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.cc#L222-L239 + + Args: + token: The S2 token string. Can be upper or lower case hex string. + + Returns: + The lat/lon (in degrees) tuple generated from the S2 token. + + Raises: + TypeError: If the token is not str. + InvalidToken: If the token length is over 16. + InvalidToken: If the token is invalid. + InvalidCellID: If the contained cell_id is invalid. + + """ + # Check input + if not token_is_valid(token): + raise InvalidToken('Cannot decode invalid S2 token: {}'.format(token)) + + # Convert to cell ID and decode to lat/lon + return cell_id_to_lat_lon(token_to_cell_id(token)) + + +# +# Token canonicalisation +# + +def token_to_canonical_token(token: str) -> str: + """ + Convert S2 token to a canonicalised S2 token. + + This produces a token that matches the form generated by the reference C++ implementation: + + - Lower case (except 'X' below) + - No whitespace + - Trailing '0' characters stripped + - Zero cell ID represented as 'X', not 'x' or '' + + Args: + token: The S2 token string to canonicalise. + + Returns: + The canonicalised S2 token. + + """ + # Convert token to lower case. + # Note that 'X' below will be returned upper case + token = token.lower() + + # Strip any surrounding whitespace + token = token.strip() + + # Strip any trailing zeros + token = token.rstrip('0') + + # If empty string or 'x', return 'X' token + if token in ('', 'x'): + token = 'X' + + return token + + +# +# Validation +# + +def cell_id_is_valid(cell_id: int) -> bool: + """ + Check that a S2 cell ID is valid. + + Looks for valid face bits and a trailing 1 bit in one of the correct locations. + + Args: + cell_id: The S2 cell integer to validate. + + Returns: + True if the cell ID is valid, False otherwise. + + Raises: + TypeError: If the cell_id is not int. + + """ + # Check input + if not isinstance(cell_id, int): + raise TypeError('Cannot decode S2 cell ID from type: {}'.format(type(cell_id))) + + # Check for zero ID + # This avoids overflow warnings below when 1 gets added to max uint64 + if cell_id == 0: + return False + + # Check face bits + if (cell_id >> _S2_POS_BITS) > 5: + return False + + # Check trailing 1 bit is in one of the even bit positions allowed for the 30 levels, using the + # mask: 0b0001010101010101010101010101010101010101010101010101010101010101 = 0x1555555555555555 + lowest_set_bit = cell_id & (~cell_id + 1) # pylint: disable=invalid-unary-operand-type + if not lowest_set_bit & 0x1555555555555555: + return False + + return True # Checks have passed, cell ID must be valid + + +def token_is_valid(token: str) -> bool: + """ + Check that a S2 token is valid. + + Looks for valid characters, then checks that the contained S2 cell ID is also valid. Note that + the '', 'x' and 'X' tokens are considered invalid, since the cell IDs they represent are + invalid. + + Args: + token: The S2 token string to validate. + + Returns: + True if the token is valid, False otherwise. + + Raises: + TypeError: If the token is not str. + + """ + # Check input + if not isinstance(token, str): + raise TypeError('Cannot check S2 token with type: {}'.format(type(token))) + + # First check string with regex + if not re.match(r'^[0-9a-fA-f]{1,16}$', token): + return False + + # Check the contained cell ID is valid + return cell_id_is_valid(token_to_cell_id(token)) + + +# +# Level extraction functions +# + +def cell_id_to_level(cell_id: int) -> int: + """ + Get the level for a S2 cell ID. + + See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.h#L543-L551 + + Args: + cell_id: The S2 cell ID integer. + + Returns: + The level of the S2 cell ID. + + Raises: + TypeError: If the cell_id is not int. + InvalidCellID: If the cell_id is invalid. + + """ + # Check input + if not cell_id_is_valid(cell_id): + raise InvalidCellID('Cannot decode invalid S2 cell ID: {}'.format(cell_id)) + + # Find the position of the lowest set one bit, which will be the trailing one bit. The level is + # given by the max level (30) minus the floored division by two of the position of the lowest + # set bit. + # + # The position of the lowest set bit is found using 'count trailing zeros', which would be + # equivalent to the C++20 function std::countr_zero() or the ctz instruction. + lsb_pos = 0 + while cell_id != 0: # pragma: no branch + if cell_id & 1: + break + lsb_pos += 1 + cell_id >>= 1 + + return int(_S2_MAX_LEVEL - (lsb_pos >> 1)) + + +def token_to_level(token: str) -> int: + """ + Get the level for a S2 token. + + See s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.h#L543-L551 + + Args: + token: The S2 token string. Can be upper or lower case hex string. + + Returns: + The level of the S2 token. + + Raises: + TypeError: If the token is not str. + InvalidToken: If the token length is over 16. + InvalidToken: If the token is invalid. + InvalidCellID: If the contained cell_id is invalid. + + """ + # Check input + if not token_is_valid(token): + raise InvalidToken('Cannot decode invalid S2 token: {}'.format(token)) + + # Convert to cell ID and get the level for that + return cell_id_to_level(token_to_cell_id(token)) + + +# +# Parent functions +# + +def cell_id_to_parent_cell_id( + cell_id: int, level: Optional[int] = None +) -> int: + """ + Get the parent cell ID of a S2 cell ID. + + Args: + cell_id: The S2 cell ID integer. + level: The parent level to get the cell ID for. Must be less than or equal to the current + level of the provided cell ID. If unspecified, or None, the direct parent cell ID will + be returned. + + Returns: + The parent cell ID at the specified level. + + Raises: + TypeError: If the cell_id is not int. + InvalidCellID: If the cell_id is invalid. + ValueError: If cell ID is already level 0 and level is None. + ValueError: When level is not an integer, is < 0 or is > 30. + ValueError: If level is greater than the provided cell ID level. + + """ + # Check input + if not cell_id_is_valid(cell_id): + raise InvalidCellID('Cannot decode invalid S2 cell ID: {}'.format(cell_id)) + + # Get current level of the cell ID and check it is suitable with the requested level + current_level = cell_id_to_level(cell_id) + if level is None and current_level == 0: + raise ValueError('Cannot get parent cell ID of a level 0 cell ID') + if level is None: + level = current_level - 1 + + if not isinstance(level, int) or level < 0 or level > _S2_MAX_LEVEL: + raise ValueError('S2 level must be integer >= 0 and <= 30') + + if level > current_level: + raise ValueError('Cannot get level {} parent cell ID of cell ID with level {}'.format( + level, current_level + )) + if level == current_level: + # Requested parent level is current level, return cell ID itself + return cell_id + + # Truncate to desired level + # This is done by finding the mask of the trailing 1 bit for the specified level, then zeroing + # out all bits less significant than this, then finally setting the trailing 1 bit. This is + # still necessary to do even after a reduced number of steps `required_steps` above, since each + # step contains multiple levels that may need partial overwrite. Additionally, we need to add + # the trailing 1 bit, which is not yet set above. + least_significant_bit_mask = 1 << (2 * (_S2_MAX_LEVEL - level)) + cell_id = (cell_id & -least_significant_bit_mask) | least_significant_bit_mask + + return cell_id + + +def token_to_parent_token(token: str, level: Optional[int] = None) -> str: + """ + Get the parent token of a S2 token. + + Args: + token: The S2 token string. Can be upper or lower case hex string. + level: The parent level to get the token for. Must be less than or equal to the current + level of the provided toke. If unspecified, or None, the direct parent token will be + returned. + + Returns: + The parent token at the specified level. + + Raises: + TypeError: If the token is not str. + InvalidToken: If the token length is over 16. + InvalidToken: If the token is invalid. + InvalidCellID: If the contained cell_id is invalid. + ValueError: If token is already level 0 and level is None. + ValueError: When level is not an integer, is < 0 or is > 30. + ValueError: If level is greater than the provided token level. + + """ + # Check input + if not token_is_valid(token): + raise InvalidToken('Cannot decode invalid S2 token: {}'.format(token)) + + # Convert to cell ID and get parent and convert back to token + return cell_id_to_token(cell_id_to_parent_cell_id(token_to_cell_id(token), level)) diff --git a/setup.py b/setup.py index 8281bd0..886966f 100644 --- a/setup.py +++ b/setup.py @@ -32,16 +32,17 @@ install_requires=[], extras_require={ 'dev': [ - 'flake8~=3.8', - 'nox~=2020.8.22', - 'pydocstyle~=5.1', - 'pylint~=2.6.0', - 'pytest~=6.2.1', - 'pytest-cov~=2.10', + 'flake8~=4.0', + 'furo==2021.11.23', + 'nox~=2021.10.1', + 'pydocstyle~=6.1.1', + 'pylint~=2.12.1', + 'pytest~=6.2.5', + 'pytest-cov~=3.0.0', 'pytest-instafail~=0.4.2', 'pytest-xdist~=2.2.0', - 'Sphinx~=3.4.0', - 'sphinx-redactor-theme==0.0.1', + 'Sphinx~=4.3.1', + 'sphinx-notfound-page~=0.8', 'sphinx-sitemap==2.2.0', ], }, diff --git a/tests/test_s2cell.py b/tests/test_s2cell.py index ea5acfc..6b9bf6e 100644 --- a/tests/test_s2cell.py +++ b/tests/test_s2cell.py @@ -20,11 +20,12 @@ import pytest import s2cell +from s2cell.s2cell import _S2_POS_BITS, _s2_face_uv_to_xyz def test_invalid__s2_face_uv_to_xyz(): with pytest.raises(ValueError, match=re.escape('Cannot convert UV to XYZ with invalid face: 6')): - s2cell._s2_face_uv_to_xyz(6, (0, 0)) + _s2_face_uv_to_xyz(6, (0, 0)) def test_zero_cell_id_to_token(): @@ -114,7 +115,7 @@ def test_invalid_cell_id_to_lat_lon(): s2cell.cell_id_to_lat_lon(1.0) with pytest.raises(s2cell.InvalidCellID, match=re.escape('Cannot decode invalid S2 cell ID: 13835058055282163712')): - s2cell.cell_id_to_lat_lon(int(0b110 << s2cell._S2_POS_BITS)) + s2cell.cell_id_to_lat_lon(int(0b110 << _S2_POS_BITS)) def test_cell_id_to_lat_lon_compat(): @@ -134,7 +135,7 @@ def test_invalid_token_to_lat_lon(): s2cell.token_to_lat_lon('a' * 17) with pytest.raises(s2cell.InvalidToken, match=re.escape('Cannot decode invalid S2 token: c000000000000000')): - s2cell.token_to_lat_lon('{:016x}'.format(0b110 << s2cell._S2_POS_BITS)) + s2cell.token_to_lat_lon('{:016x}'.format(0b110 << _S2_POS_BITS)) def test_token_to_lat_lon_compat(): @@ -166,9 +167,9 @@ def test_token_to_canonical_token(token, expected): (0b1111010101010101010101010101010101010101010101010101010101010101, False), # Invalid face (0, False), ] + [ - (1 << even_number, True) for even_number in range(0, s2cell._S2_POS_BITS, 2) + (1 << even_number, True) for even_number in range(0, _S2_POS_BITS, 2) ] + [ - (1 << odd_number, False) for odd_number in range(1, s2cell._S2_POS_BITS, 2) + (1 << odd_number, False) for odd_number in range(1, _S2_POS_BITS, 2) ]) def test_cell_id_is_valid(cell_id, is_valid): assert s2cell.cell_id_is_valid(cell_id) == is_valid