Created
July 29, 2020 13:10
-
-
Save CloudNiner/e5552aaae6c235193b8f196558b92e2c to your computer and use it in GitHub Desktop.
Query OSMX by S2Region in Python
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import logging | |
import time | |
import osmx | |
import pywraps2 as s2 | |
logging.basicConfig(format="%(message)s") | |
logger = logging.getLogger(__name__) | |
logger.setLevel(logging.DEBUG) | |
OSMX_S2_CELL_MAX_CELLS = 1024 | |
OSMX_S2_CELL_MAX_LEVEL = 16 | |
def traverse_cell(cell_id, cell_level): | |
""" Traverse all child cells for cell_id at the desired cell_level. """ | |
start = cell_id.child_begin(cell_level) | |
end = cell_id.child_end(cell_level) | |
current_cell = start | |
while current_cell.id() != end.id(): | |
yield current_cell | |
current_cell = current_cell.next() | |
def node_ids_for_region(osmx_txn, s2_region, min_cell_level=None): | |
""" Get node_ids from osmx within the s2_region. """ | |
start_time = time.time() | |
cell_node_index = osmx.Index(osmx_txn, b"cell_node") | |
coverer = s2.S2RegionCoverer() | |
# Values from https://github.com/protomaps/OSMExpress/blob/master/src/extract.cpp#L131-L132 | |
coverer.set_max_cells(OSMX_S2_CELL_MAX_CELLS) | |
if isinstance(min_cell_level, int): | |
coverer.set_min_level(min_cell_level) | |
coverer.set_max_level(OSMX_S2_CELL_MAX_LEVEL) | |
s2_cell_ids = coverer.GetCovering(s2_region) | |
logger.debug("Region covers {} cells.".format(len(s2_cell_ids))) | |
node_ids = list() | |
for cell_id in s2_cell_ids: | |
for max_level_cell_id in traverse_cell(cell_id, OSMX_S2_CELL_MAX_LEVEL): | |
node_ids.extend(cell_node_index.get(max_level_cell_id.id())) | |
node_id_set = set(node_ids) | |
logger.debug( | |
"Returned {} node_ids in {}s".format(len(node_id_set), time.time() - start_time) | |
) | |
return node_id_set |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
""" Convert various geometries to S2 regions | |
This method isn't used in this example gist, but it could be used to convert a polygon | |
for input to search_by_tags, e.g.: | |
osmx_db = "/path/to/db.osmx" | |
polygon = ??? | |
tags = ["name=", "highway=primary"] | |
features = search_by_tags(osmx_db, polygon_to_s2_region(polygon), tags) | |
""" | |
import pywraps2 as s2 | |
def polygon_to_s2_region(polygon): | |
""" Convert a geojson.Polygon to s2.S2Region. """ | |
loops = [ | |
s2.S2Loop( | |
list( | |
map( | |
lambda coord: s2.S2LatLng.FromDegrees(coord[1], coord[0]) | |
.Normalized() | |
.ToPoint(), | |
loop[:-1], | |
) | |
) | |
) | |
for loop in polygon.coordinates | |
] | |
for loop in loops: | |
loop.Normalize() | |
return s2.S2Polygon(*loops) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment