Returns an s2 covering as ndgeojson

s2cover.py

#!/usr/bin/env python3
import json
import argparse
import s2sphere as s2


def point2coord(point):
    return s2.LatLng.longitude(point).degrees, s2.LatLng.latitude(point).degrees


def cellid2feature(cellid):
    cell = s2.Cell(cellid)
    coords = [point2coord(cell.get_vertex(i)) for i in range(4)]
    return {
        "type": "Feature",
        "properties": {
            "id": cellid.id()
        },
        "geometry": {
            "type": "Polygon",
            "coordinates": [coords + coords[:1]]
        }
    }


def bbox2region(bbox):
    a = s2.LatLng.from_degrees(bbox[1], bbox[0])
    b = s2.LatLng.from_degrees(bbox[3], bbox[2])
    return s2.LatLngRect.from_point_pair(a, b)


def main():
    parser = argparse.ArgumentParser(description="Return s2 covering as ndgeojson")
    parser.add_argument('bbox', type=float, nargs=4)
    parser.add_argument('-m', '--min-level', default=14, type=int)
    parser.add_argument('-M', '--max-level', default=14, type=int)
    args = parser.parse_args()

    region = bbox2region(args.bbox)
    rc = s2.RegionCoverer()
    rc.max_level, rc.min_level = args.max_level, args.min_level

    for cid in rc.get_covering(region):
        print(json.dumps(cellid2feature(cid)))


if __name__ == '__main__':
    main()