Creating shapefile of S2 cells for given level

Question

I want to create a shapefile that has polygons matching with the outlines of the S2 cells (http://s2geometry.io/) of a given level, within a certain bounding box.

Roughly what this website does: Sidewalk Labs Region Coverer, but giving me the data as an actual shapefile.

I'll accept any working option:

• generated through a website

• though QGIS (and/or GRASS) if properly explained how to (as in if you're giving me a Python script, please tell me HOW to use it within QGIS, rather then just listing me all the options for the parameters)

Whatever works and is understandable and doable without having to purchase the software.

---

EDIT: if it requires hard-coding the S2 cell level, I'm particularly interested in levels 14 and 17, if the code requires manual altering the level please make sure it's clearly visible what value needs editing.

Answer 1

I came across this thread the other day with the identical need, and after searching around a bunch, couldn't find anything ready-made. After biting the bullet and sitting down and playing a bit with the Python S2Sphere library a bit, it turned out to be more straightforward than I'd hoped, and I banged out the following bit of code in a couple of hours.

It's a Python 3 program. (s2sphere doesn't want to import properly in Python 2, at least with my machine's setup. This means that if you are using QGIS 2.x, you cannot use QGIS' bundled Python to run it, you will need a separate Python 3 installation).

(You will also, of course, need to install the aforementioned s2sphere library into your Python 3 installation - instructions on its homepage.)

The meat of the program is the s2gen() function. You can either import this into another python program, or you can run the whole kit and caboodle from the command line as a standalone program. But the above function and the standalone program take the same 5 arguments:

• max_lat, min_lat, max_lng, min_lng :: float (so "N, S, E, W" order)
• s2_lvl :: int

And both return a KML file containing the all s2 cells (as polygons) at the given cell level that cover the cap (i.e.: "rectangle") bounded by the given N, S, E and W lat and lng coordinates. The function returns the KML file contents as a string, and the standalone just prints them to stdout (so you can redirect into a KML file of your choosing). Then this KML can of course be imported into QGIS or Google Earth or whatever.

So for example (if you save the code as s2gen.py) if you wanted to generate a KML containing the level 14 S2 cell polygons over New York's Central Park, then first you'd find the N/S/E/W extents, which happen to be

• max_lat (N): 40.800457
• min_lat (S): 40.764412
• max_lng (E): -73.949232
• min_lng (W): -73.981716

And then just dump the coordinates into the script, redirecting the output to your file:

C:\...\> python3 s2gen.py 40.800457 40.764412 -73.949232 -73.981716 14 > cpark_s2l14.kml

The KML file being one that you can import into QGIS.

Or within python:

from s2gen import s2gen

kml = s2gen(40.800457, 40.764412, -73.949232, -73.981716, 14)

A caution that because the number of cells grows exponentially per s2 level, the files generated can get pretty big, depending on level and area. I was generating grids to cover all of the city of Calgary, which is about 825 km², and I pretty much maxed out at level 19. Anything more than that was too big for QGIS to import without crashing.

Anyway, the code. Again, save this as s2gen.py:

# -*- coding: utf-8 -*-

# Simple s2 KML generator
# Copyright 2018, Sean C. Nichols
#  [email protected]

# Hereby released into the public domain: feel free to abuse this code as you see fit!
#

# usage: python s2gen.py <max_lat> <min_lat> <max_lng> <min_lng> <s2_lvl>
#
# (i.e.: N-S-E-W order)
#

from __future__ import print_function

from sys import argv, stderr

from s2sphere import RegionCoverer, Cell, LatLng, LatLngRect

def usage(msg=None):
    if msg:
        stderr.write(msg + '\n\n')
    stderr.write('usage: {} <max_lat> <min_lat> <max_lng> <min_lng> <s2_lvl>\n'.format(argv[0]))
    exit(1)

def swap_latlong(latlong):
    return '{},{}'.format(*latlong.split(',')[::-1])

def s2gen(max_lat, min_lat, max_lng, min_lng, s2_lvl):
    point_nw = LatLng.from_degrees(max_lat, min_lng)
    point_se = LatLng.from_degrees(min_lat, max_lng)

    rc = RegionCoverer()
    rc.min_level = s2_lvl
    rc.max_level = s2_lvl
    rc.max_cells = 1000000

    cellids = rc.get_covering(LatLngRect.from_point_pair(point_nw, point_se))

    kml = (
        '<?xml version="1.0" encoding="UTF-8"?>\n'
        '<kml xmlns="http://www.opengis.net/kml/2.2" xmlns:kml="http://www.opengis.net/kml/2.2" xmlns:atom="http://www.w3.org/2005/Atom">\n'
        '<Document><name>{0}</name>\n'
        '<Style id="s2_poly_style"><LineStyle><color>ff0000ff</color><width>2</width></LineStyle><PolyStyle><fill>0</fill></PolyStyle></Style>\n'
        '<Folder><name>{0}</name><open>1</open>'
    ).format('S2 cells level ' + str(s2_lvl))
    for cid in cellids:
        vertices = [LatLng.from_point(Cell(cid).get_vertex(v)) for v in range(4)]
        kml_coords = ['{},0'.format(swap_latlong(str(v).split()[-1])) for v in vertices]
        kml += (
            '<Placemark><name>{}</name><styleUrl>#s2_poly_style</styleUrl><Polygon><tessellate>1</tessellate><outerBoundaryIs><LinearRing>'
            '<coordinates>{} {}</coordinates>'
            '</LinearRing></outerBoundaryIs></Polygon></Placemark>'
        ).format(cid.id(), ' '.join(kml_coords), kml_coords[0])
    kml += (
        '</Folder></Document>\n'
        '</kml>'
    )

    return kml

if __name__ == '__main__':
    if len(argv) != 6:
        usage()

    try:
        (max_lat, min_lat, max_lng, min_lng) = map(float, argv[1:5])
        s2_lvl = int(argv[5])
    except ValueError:
        usage('Arguments must be numbers.')

    print(s2gen(max_lat, min_lat, max_lng, min_lng, s2_lvl))

Answer 2

@mrputter's solution is a good start, but I want to use the s2cells as a grid for snapping purposes so I don't need the full coordinates for every cell. With that in mind, I can afford to just store the first vertex (e.g. bottom left corner), so instead of storing 4 lines (=4 pairs of points) each cell can be represented by one point. Hence, his solution can be adapted to work for cells up to at least lvl 22 (since the output size is reduced by factor of 8 (=2^3)).

def s2cover(max_lat, min_lat, max_lng, min_lng, s2_lvl):
    point_nw = LatLng.from_degrees(max_lat, min_lng)
    point_se = LatLng.from_degrees(min_lat, max_lng)

    rc = RegionCoverer()
    rc.min_level = s2_lvl
    rc.max_level = s2_lvl
    rc.max_cells = 1000000

    cellids = rc.get_covering(LatLngRect.from_point_pair(point_nw, point_se))
    return cellids

def s2geojson(cellids):
    fc = {
      "type": "FeatureCollection"
    }

    features = []

    for cid in cellids:
        pt_feature = {
          "type": "Feature",
          "properties": {},
          "geometry": {
            "type": "Point"
          }
        }
        cell = Cell(cid)
        latlong = LatLng.from_point(cell.get_vertex(0))
        pt_feature['geometry']['coordinates'] = [float(s) for s in str(latlong).split()[-1].split(',')[::-1]]
        pt_feature['properties']['s2_cellid'] = cid.id()

        features.append(pt_feature)

    fc['features'] = features

    return fc

As an added bonus, there is minimal distortion at the higher zoom levels, so if there is a need to draw the grid outlines, we can use the points at the perimeter to draw a minimal number of lines (total number of rows + total number of columns).