BabelGrid is a common python API to work with different established geospatial indexing systems.
Currently, it supports H3, S2 and Bing geospatial indexing systems. BabelGrid does not have the intention to replace any of the existing APIs, but to create a common framework that geo-folks can use to easily switch between grids. Also, it generates methods around the tiles that ease the data analysis pipeline with seamlessly integration with well knonw libraries such as Shapely and GeoPandas.
Install with
pip install babelgrid
>>> from babelgrid import Babel
tile = Babel('h3').geo_to_tile(lat=-23, lon=-43, area_km=1) tile Tile: grid_type "h3", resolution 8, tile_id 88a8a2b66dfffff
Access the geojson, wkt and shapely descriptions of the tile:
>>> tile.geometry.geojson
{'type': 'Polygon',
'coordinates': (((-42.99741709893686, -23.004282833594505),
(-42.9932470321478, -23.00127887552568),
(-42.994161748920796, -22.996608473771282),
(-42.99924646130203, -22.994942061847414),
(-43.00341650043048, -22.997946087213307),
(-43.002501854850166, -23.002616457194414),
(-42.99741709893686, -23.004282833594505)),)}
>>> tile.geometry.wkt
'POLYGON ((-42.9974170989368574 -23.0042828335945053, -42.9932470321477993 -23.0012788755256814, -42.9941617489207957 -22.9966084737712819, -42.9992464613020289 -22.9949420618474143, -43.0034165004304825 -22.9979460872133075, -43.0025018548501663 -23.0026164571944136, -42.9974170989368574 -23.0042828335945053))'
>>> tile.geometry.shapely
Fill a geometry with s2 tiles of resolution 10
>>> tiles = Babel('s2').polyfill(geometry, resolution=10)
>>> tiles
[Tile: grid_type "s2", resolution 10, tile_id 94d28d,... ,Tile: grid_type "s2", resolution 10, tile_id 94d28f]
Load a geopandas dataframe with the selected tiles
>>> import geopandas as gpd
>>> gpd.GeoDataFrame([t.to_dict() for t in tiles], geometry='shapely')
You have to initialize the Babel object with any of the available grids.
>>> Babel.available_grids()
['s2', 'h3', 'bing']
>>> grid = Babel('s2') # example
It receives a coordinate pair (lat, lon) and either the native grid resolution or an area in km2. If it receives an area, it automatically finds what is the resolution for that tile system and latitute that best approximates the given area.
>>> Babel('s2').geo_to_tile(2, 3, resolution=10)
Tile: grid_type "s2", resolution 10, tile_id 100fb1
>>> Babel('bing').geo_to_tile(2, 3, area_km=0.1)
Tile: grid_type "bing", resolution 17, tile_id 12222230201200322
>>> Babel('bing').geo_to_tile(2, 3, area_km=0.1).area_km
0.0934819087
It receives a tile id and converts it to a Tile Object.
>>> Babel('s2').id_to_tile('100fb1')
Tile: grid_type "s2", resolution 10, tile_id 100fb1
One of the most common uses to geospatial indexing systems is to fill up a geometry. This function receives a geometry that can be a polygon or multipolygons and returns a list of Tile Objects.
>>> tiles = Babel('s2').polyfill(geometry, resolution=10)
>>> tiles
[Tile: grid_type "s2", resolution 10, tile_id 94d28d,... ,Tile: grid_type "s2", resolution 10, tile_id 94d28f]
You can also pass a 'desired' grid area using the parameter grid_km
.
>>> tiles = Babel('bing').polyfill(geometry, area_km=10)
>>> tiles
[Tile: grid_type "bing", resolution 14, tile_id 21031113121331, ..., Tile: grid_type "bing", resolution 14, tile_id 21031113121333]
The image below shows polyfill
being applied for the same geometry for different grid types and sizes.
The Tile Object is a central piece of the package. This is the object that is returned by most of the methods implemented. It is good because it has some handy features that speed-up the analysis process.
>>> tile.geometry.wkt
>>> tile.geometry.geojson
>>> tile.geometry.shapely
>>> tile.to_parent()
>>> tile.to_children()
>>> tile.area_km
>>> tile.to_dict()
H3 | S2 | BING/QuadTree | |
---|---|---|---|
Tile Shape | Hexagonal | Square | Square |
Resolution Range | 0 - 15 | 0 - 30 | 1 - 23 (infinite) |
API Reference | h3-py | s2sphere | pygeotile |
Original Documentation | H3 | S2 Geometry | Bing Maps Tile System |
:star: Kudos to all developer of H3, S2 and Bing/QuadTree systems.
Lookup table with grid resolutions at equator by area in km2. Note that the area is written in scientific notation (10^x) and x is the index of the table.
Area (10^x km2) | H3 | S2 | BING/QuadTree |
---|---|---|---|
9 | - | - | 1 |
8 | - | 0 | 2 |
7 | - | 1,2 | 3,4 |
6 | 0,1 | 3,4 | 5,6 |
5 | 2 | 5 | 7 |
4 | 3 | 6,7 | 8,9 |
3 | 4 | 8 | 10,11 |
2 | 5 | 9,10 | 12 |
1 | 6,7 | 11,12 | 13,14 |
0 | 8 | 13 | 15,16 |
-1 | 9 | 14,15 | 17 |
-2 | 10 | 16,17 | 18,19 |
-3 | 11 | 18 | 20,21 |
-4 | 12,13 | 19,20 | 22 |
-5 | 14 | 21,22 | 23 |
-6 | 15 | 23 | - |
-7 | - | 24,25 | - |
-8 | - | 26,27 | - |
-9 | - | 28 | - |
-10 | - | 29,30 | - |
Depending on how the tile system is built, the area of the tile varies given the latitude. For inter-region comparissons, this behaviour can affect the analysis.
The figure below shows the tile area distortion by geospatial indexing systems. The distortion is defined as
where is the tile area and the area given a latitude and the equator area. The figure shows the mean distortion given all resolutions and the error bar is the standard deviation.
Any contribution is very welcomed. You can contribute in several ways:
Start envorinment with
make create-env
Update envorinment with
make update-env
Publish to PyPi
poetry version [patch, minor, major]
make publish
This work is licensed under AM-331-A3 - see the LICENSE.md file for details.