Google S2 with Python & Jupyter

Jan 18, 2018 - Geo Python

Google is working again on S2 a spatial library !!!

And they even have created a website to communicate about it s2geometry.

The C++ port contains a Python Swig interface.

I’ve been using an unofficial Python port with Jupyter for years now things are way more simpler.

If you are on Arch I’ve create a package, simply install AUR s2geometry-git

First we want a clean Jupyter install from scratch:

virtualenv3 ~/dev/venv3
source ~/dev/venv3/bin/activate
pip install jupyter
pip install cython
pip install numpy
pip install matplotlib scikit-learn scipy Shapely folium geojson Cartopy
cp /usr/lib/python3.6/site-packages/_pywraps2.so $VIRTUAL_ENV/lib/python3.6/site-packages      
cp /usr/lib/python3.6/site-packages/pywraps2.py $VIRTUAL_ENV/lib/python3.6/site-packages

Here is a simple test map.

import folium
import pywraps2 as s2

# create a rect in s2
region_rect = s2.S2LatLngRect(
        s2.S2LatLng.FromDegrees(48.831776, 2.222639),
        s2.S2LatLng.FromDegrees(48.902839, 2.406))

# ask s2 to create a cover of this rect
coverer = s2.S2RegionCoverer()
coverer.set_min_level(10)
coverer.set_max_level(30)
coverer.set_max_cells(60)
covering = coverer.GetCovering(region_rect)
print([c.ToToken() for c in covering])

# create a map
map_osm = folium.Map(location=[48.86, 2.3],zoom_start=12, tiles='Stamen Toner')

# get vertices from rect to draw them on map
rect_vertices = []
for i in [0, 1, 2, 3, 0]:
    vertex = region_rect.GetVertex(i)
    rect_vertices.append([vertex.lat().degrees(), vertex.lng().degrees()])
    
# draw the cells
style_function = lambda x: {'weight': 1, 'fillColor':'#eea500'}
for cellid in covering:
    cell = s2.S2Cell(cellid)
    vertices = []
    for i in range(0, 4):
        vertex = cell.GetVertex(i)
        latlng = s2.S2LatLng(vertex)
        vertices.append([latlng.lng().degrees(),
                         latlng.lat().degrees()])
        
    gj = folium.GeoJson({ "type": "Polygon", "coordinates": [vertices]}, style_function=style_function)
    gj.add_children(folium.Popup(cellid.ToToken()))
    gj.add_to(map_osm)
    
# warning PolyLine is lat,lng based while GeoJSON is not
ls = folium.PolyLine(rect_vertices, color='red', weight=2)
ls.add_children(folium.Popup("shape"))
ls.add_to(map_osm)
    
map_osm

And here is the resulting Jupyter notebook

See Also