aaliddell / s2cell Goto Github PK

View Code? Open in Web Editor NEW

27.0 2.0 3.0 11.75 MB

Minimal Python S2 Geometry cell ID, token and lat/lon conversion library

Home Page: https://docs.s2cell.aliddell.com

License: Apache License 2.0

Python 69.69% Makefile 0.07% C++ 4.27% HTML 25.97%

s2 python python3 geometry spherical-geometry

s2cell's Introduction

s2cell

Minimal Python S2 Geometry cell ID, token and lat/lon conversion library.

Docs | PyPI | GitHub

Overview

This library does conversion between S2 cell ID, S2 token and latitude/longitude and was written as a method to understand the way the S2 cell system works; hopefully this is useful to others as a single-file reference on the process, where tracing the relevant parts from the reference C++ implementation can be somewhat daunting. All steps in the conversions are well commented and written to be understandable and functional rather than strictly fast, although little is different from the reference implementation.

The library is checked against a test suite generated from the reference C++ implementation to ensure conformity with the standard.

Should you need more complete S2 Geometry functionality or a fast C-based implementation, please consider using the Python extension included in the s2geometry repository or the pure-Python s2sphere implementation.

Issues and PRs are very welcome to improve the implementation, descriptions or to correct any misunderstandings of the S2 cell system. Please note that this library strives to be an easy to read reference rather than aiming for peak performance (it is in Python after all), so PRs which reduce readability of the implementation (such as for Python specific speed optimisations) are generally discouraged. However, if you have optimisations that are applicable to S2 implementations across many languages and can be described easily, please do consider making a PR.

Installation

This package can be installed from PyPI with pip or any other Python package manager:

pip install s2cell

Usage

The full documentation, including the API Reference, is available at docs.s2cell.aliddell.com. Below is a quick start guide for the most common uses.

The library is designed to be minimal, predictable and purely functional. Conversion from lat/lon (in degrees) to a cell ID or token can be done with the following two functions:

s2cell.lat_lon_to_cell_id(-10.490091, 105.641318)  # -> 3383782026967071427
s2cell.lat_lon_to_token(-10.490091, 105.641318)    # -> '2ef59bd352b93ac3'

By default, these conversions will give you a level 30 leaf-cell as output. If you require a lower precision level, you can specify this:

s2cell.lat_lon_to_cell_id(-10.490091, 105.641318, 10)  # -> 3383781119341101056
s2cell.lat_lon_to_token(-10.490091, 105.641318, 0)     # -> '3'

Conversion from a cell ID or token to lat/lon (in degrees) can be done with the following two functions:

s2cell.cell_id_to_lat_lon(3383781119341101056)  # -> (-10.452552407574101, 105.6412526632361)
s2cell.token_to_lat_lon('3')                    # -> (0.0, 90.0)

The lat/lon returned will be the center of the cell at the level available in the provided cell ID or token.

There are also a few other useful functions for inspecting or converting a cell ID/token:

# Conversion between cell ID and token
s2cell.cell_id_to_token(3383781119341101056)  # -> '2ef59b'
s2cell.token_to_cell_id('3')                  # -> 3458764513820540928

# Level extraction
s2cell.cell_id_to_level(3383781119341101056)  # -> 10
s2cell.token_to_level('3')                    # -> 0

# Parent cell calculation
s2cell.cell_id_to_parent_cell_id(3383781119341101056)     # -> 3383782218852728832
s2cell.cell_id_to_parent_cell_id(3383781119341101056, 2)  # -> 3386706919782612992

s2cell.token_to_parent_token('2ef59b')                    # -> '2ef59c'
s2cell.token_to_parent_token('2ef59b', 2)                 # -> '2f'

# Token canonicalisation
s2cell.token_to_canonical_token('2ef59BD352b90') # -> '2ef59bd352b9'

Useful S2 Geometry Links

A list of useful links for S2 related concepts and projects can be found here.

License

This project is released under the same license as the reference C++ S2 Geometry implementation, namely the Apache 2.0 License.

s2cell's People

Contributors

Stargazers

Watchers

Forkers

tdunning geopipe igorgatis

s2cell's Issues

Fix docs line spacing

Paragraphs and sections have insufficient spacing. Perhaps consider switching themes...

Add alternative to http://s2.sidewalklabs.com/regioncoverer/

The http://s2.sidewalklabs.com/regioncoverer/ has been down for a while now. You may want to add the following links to the list of visualizations:

https://s2.inair.space/ whose related post worth mentioning https://blog.nobugware.com/post/2018/wasm_go_s2_javascript/
https://igorgatis.github.io/ws2/ -my fork from the link above which has some of the functionality that sidewalklabs' provided.

Add support for neighborhood generation

It would be really nice if there was a function to generate neighborhood cells for a given cell id. More specifically:

I managed to hack these functions reusing some of the current code. It seems that it would require a bit of refactor to make function reusable. I attached my code down below (which I can send a PR if you are willing to accept it).

I believe there is something wrong for the cube corner case. I was expecting a tile right here link:

import math

from s2cell import s2cell  # type: ignore
from s2cell.s2cell import (  # type: ignore
    _S2_INVERT_MASK,
    _S2_LOOKUP_BITS,
    _S2_MAX_LEVEL,
    _S2_MAX_SIZE,
    _S2_POS_BITS,
    _S2_SWAP_MASK,
    _s2_face_uv_to_xyz,
    _s2_st_to_ij,
    cell_id_to_level,
)


def get_edge_neighbor_cell_ids(cellid: int) -> list[int]:
    level = cell_id_to_level(cellid)
    size = _get_size_ij(level)
    face, i, j = _to_face_ij_orientation(cellid)
    return [
        _from_face_ij_same(face, i, j - size, j - size >= 0, level),
        _from_face_ij_same(face, i + size, j, i + size < _S2_MAX_SIZE, level),
        _from_face_ij_same(face, i, j + size, j + size < _S2_MAX_SIZE, level),
        _from_face_ij_same(face, i - size, j, i - size >= 0, level),
    ]


def get_vertex_neighbor_cell_ids(cellid: int) -> list[int]:
    level = cell_id_to_level(cellid)
    size = _get_size_ij(level)
    face, i, j = _to_face_ij_orientation(cellid)
    neighbors = []
    for di in (-size, size):
        for dj in (-size, size):
            same = (0 <= (i + di) < _S2_MAX_SIZE) and (0 <= (j + dj) < _S2_MAX_SIZE)
            neighbors.append(_from_face_ij_same(face, i + di, j + dj, same, level))
    return neighbors


def _get_size_ij(level: int) -> int:
    return 1 << (_S2_MAX_LEVEL - level)


def _to_face_ij_orientation(cell_id: int) -> tuple[int, int, int]:
    s2cell._s2_init_lookups()

    # https://github.com/aaliddell/s2cell/blob/main/s2cell/s2cell.py#L546-L607
    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 = s2cell._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
    return face, i, j


def _get_face(p: tuple[float, float, float]) -> int:
    # https://github.com/aaliddell/s2cell/blob/main/s2cell/s2cell.py#L389-L391
    face = max(enumerate(p), key=lambda p: abs(p[1]))[0]  # Largest absolute component
    if p[face] < 0.0:
        face += 3
    return face


def _xyz_to_face_uv(p: tuple[float, float, float]) -> tuple[int, float, float]:
    face = _get_face(p)
    # https://github.com/aaliddell/s2cell/blob/main/s2cell/s2cell.py#L393-L423
    uv = (  # pylint: disable=invalid-name
        p[1 - ((face + 1) >> 1)] / p[face % 3],  # U
        p[2 - (face >> 1)] / p[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
    return face, uv[0], uv[1]
    # if face == 0:
    #     return face, p[1] / p[0], p[2] / p[0]
    # elif face == 1:
    #     return face, -p[0] / p[1], p[2] / p[1]
    # elif face == 2:
    #     return face, -p[0] / p[2], -p[1] / p[2]
    # elif face == 3:
    #     return face, p[2] / p[0], p[1] / p[0]
    # elif face == 4:
    #     return face, p[2] / p[1], -p[0] / p[1]
    # else:
    #     return face, -p[1] / p[2], -p[0] / p[2]


def _from_face_ij_same(face: int, i: int, j: int, same: bool, level: int) -> int:
    if same:
        return _from_face_ij(face, i, j, level)
    else:
        return _from_face_ij_wrap(face, i, j, level)


def _from_face_ij_wrap(face: int, i: int, j: int, level: int) -> int:
    # https://github.com/google/s2geometry/blob/c59d0ca01ae3976db7f8abdc83fcc871a3a95186/src/s2/s2cell_id.cc#L446-L477
    i = max(-1, min(_S2_MAX_SIZE, i))
    j = max(-1, max(_S2_MAX_SIZE, j))

    scale = 1.0 / _S2_MAX_SIZE
    limit = math.nextafter(1, 2)
    u = max(-limit, min(limit, scale * ((i << 1) + 1 - _S2_MAX_SIZE)))
    v = max(-limit, min(limit, scale * ((j << 1) + 1 - _S2_MAX_SIZE)))

    face, u, v = _xyz_to_face_uv(_s2_face_uv_to_xyz(face, (u, v)))
    return _from_face_ij(
        face, _s2_st_to_ij(0.5 * (u + 1)), _s2_st_to_ij(0.5 * (v + 1)), level
    )


def _from_face_ij(face: int, i: int, j: int, level: int) -> int:
    s2cell._s2_init_lookups()

    # https://github.com/aaliddell/s2cell/blob/main/s2cell/s2cell.py#L433-L484
    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 += ((i >> offset) & lookup_mask) << (_S2_LOOKUP_BITS + 2)
        bits += ((j >> offset) & lookup_mask) << 2

        # Map bits from iiiijjjjoo to ppppppppoo using lookup table
        bits = s2cell._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

if __name__ == "__main__":
    # lat, lng = 0, 0
    # lat, lng = -7.8556806, -34.8624166
    lat, lng = 36.7066928, -45.0000000
    level = 5
    from s2cell import lat_lon_to_cell_id, cell_id_to_token
    cell_id = lat_lon_to_cell_id(lat, lng, level)
    print(cell_id_to_token(cell_id))
    for neighbor in get_edge_neighbor_cell_ids(cell_id):
        print("  ", cell_id_to_token(neighbor))
    for neighbor in get_vertex_neighbor_cell_ids(cell_id):
        print("  ", cell_id_to_token(neighbor))

s2 in kotlin

https://github.com/Enovea/s2-geometry-kotlin

How to find neighbors

Hi. Thank you for your work in this library. I have a feature request:

Given a reference point, and a heterogeneous set of points, one might want to search some points in the set based on the distance with the reference point. In the geohash world this is usually done by computing the 8 neighbor cells.

A ticket at s2geometry suggest to build a S2PointIndex and querying it with a S2ClosestPointQuery (in the cpp implementation). It seems those datastructures are not present in this library, and this is fine, since I would prefer to use my own datastructures (or at least the one my database library allow me to use).

I suggest implementing a method that, given a cell and a precision, returns all the neighbor cells.

This would allow me to implement the algorithm described in this SO answer.

It is possible that there are some concepts about S2 that I do not really get though. In this case maybe a word about this in the documentation would be welcome?

What do you think?

Better tokens suggestion

Hi @aaliddell, good work! s2cell is a good project to test new efficient geocodes. This issue is a concept-proof suggestion.

As you say at s2-tokens,

S2 tokens can be considered analogous to the Geohash encoding system, (...). However, unlike Geohash, you cannot just chop off characters from a high precision S2 token to get a parent lower precision token, since the trailing 1 bit in the cell ID would not be set correctly in most cases.

Well, it's a "one bit problem", let's fix this bit! All the cell ID, except face prefix, is a good hiearchical code: each 2 bits (so base4 number) represents a hierarchical level... I think we can fix it by answering "Do we need all 32 bits for all applications?"

There are some application for the level30 bit? it is a cell of ~8 mm of side (0.008 m!)... Let's check a real life application, as addresses: the uncertainty ranges from 3 meters (urban areas) to 15 meters (rural areas), as illustrated below.

... So we can discard a bit! We can do a bitwise operation of rotate right in the 32 bits integer representation, interpretating it as a adding a 0-bit to the face, resulting in 4 bits instead 3 bits as face identifier.

Using your level-1 example, 00101100...000, after rotate will be 00010110...000. Your level29 example, 001011101...00100, will be 0001011101...0010.

We can encode/decode this cell ID representation (the "hierarchical token" - htk) to original token (otk) using rotate:

otk2htk(oID) rotate right
htk2otk(hID) rotate left

What do you think?

Move docs to readthedocs

Current solution does not allow browsing old docs versions and is a hack on top of GH pages

Dealing with invalid cellid

Hi all,

Without getting into details, the way I'm using the cell_id causes them to be slightly altered.

Because of this, they fail the check at cell_id_is_valid(cell_id: int) -> bool:

Is there a simple way to convert an invalid cell_id into a nearby valid cell_id?

An example of an invalid cell_id that I've created is 9926595740965789696 . (This is default level 30)

Thanks!

Add cell_id_to_vertices/token_to_vertices

Return a tuple of the coordinates of the four vertices. See https://github.com/google/s2geometry/blob/bdaaf97c60b3e29c0eb74dbdc66a7a19f1c937f6/src/s2/s2cell.h#L87-L90 and deeper

token_is_valid has typo in regex

See here

The problem is that the string "86R" passes as valid because the last character range is listed as A-f instead of A-F.