trouble generating cell ids for MultiPollygons

[kgeographer]

I'm working with a clone of the s2geometry repo as of 22 July. I was able to build a wheel for Python 3.10 on an M1 Macbook, and perform several functions with it. My goal is to generate lists of cell ids at resolution of 13 for about 17000 US counties, using MultiPolygon geometries as encoded in a Django postgresql database.

The code I have tried is below. The result for any of my county MultiPolygon objects is: ['1', '3', '5', '7', '9', 'b'], which I understand to be among the top level cells. I guess my first question is, can I expect the Python port of s2 to handle MultiPolygons, and if it can, is there anything in this code that indicates what the problem is?

Yes, this code was arrived at with the assist of ChatGPT, with many iterations regarding loop-building getting to this state.

#### ********* ####
def calculate_s2_cell_ids(polygon, max_level=13):
    # Create a region coverer
    coverer = s2.S2RegionCoverer()
    coverer.set_max_level(max_level)

    # Get the cell ids that cover the polygon
    covering = coverer.GetCovering(polygon)

    # Return the cell ids as tokens
    return [cell_id.ToToken() for cell_id in covering]

@transaction.atomic  # make the updates in a single database transaction
def update_s2_cell_ids(qs):
  for obj in qs:
    # Convert the WKT to an S2Polygon
    geom = GEOSGeometry(obj.geom)
    polygon = s2.S2Polygon()

    for part in geom:  # Iterate over each part of the multipolygon
      loop = []
      for lng, lat in part[0]:  # Outer ring
        loop.append(s2.S2LatLng.FromDegrees(lat, lng).ToPoint())
      polygon.InitNested([s2.S2Loop(loop)])  # InitNested takes a list of loops

      for inner in part[1:]:  # Inner rings (if any)
        loop = []
        for lng, lat in reversed(inner):
          loop.append(s2.S2LatLng.FromDegrees(lat, lng).ToPoint())
          # print('loop', loop)
        polygon.InitNested([s2.S2Loop(loop)])  # InitNested takes a list of loops

    print('area', polygon.GetArea())
    # Calculate the S2 cell IDs
    cell_ids = calculate_s2_cell_ids(polygon)

    print('cell_ids', cell_ids)
    # Update the s2 field (not yet)
    # obj.s2 = cell_ids
    # obj.save(update_fields=['s2'])

# 
qs = PlaceGeom.objects.filter(place__dataset='us_counties_3890')[:20]
update_s2_cell_ids(qs)

[jmr]

Your loops may be inverted. Check the areas. Are they closer to 0 or 4π?

[kgeographer]

Thank you for the quick response. For each obj, area was indeed 4π:

area 12.56632984251453
cell_ids ['1', '3', '5', '7', '9', 'b']

A Normalize() adjustment has fixed it.
Now, e.g.

area 4.07718446426448e-05
area (sqkm) 1654.9145369527257
cell_ids ['87cb9', '87cba4', '87cbbc', '87cbc4', '87cbf', '87cc04', '87cc0c', '87cc74']

the working code:

def update_s2_cell_ids(qs):
  for obj in qs:
    # Convert the WKT to an S2Polygon
    geom = GEOSGeometry(obj.geom)
    polygon = s2.S2Polygon()

    for part in geom:  # Iterate over each part of the multipolygon
      outer_loops = []
      inner_loops = []
      for lng, lat in part[0]:  # Outer ring
        outer_loops.append(s2.S2LatLng.FromDegrees(lat, lng).ToPoint())

      loop = s2.S2Loop(outer_loops)
      loop.Normalize()  # Normalize to ensure it represents the smaller region
      polygon.InitNested([loop])  # InitNested takes a list of loops

      for inner in part[1:]:  # Inner rings (if any)
        for lng, lat in reversed(inner):
          inner_loops.append(s2.S2LatLng.FromDegrees(lat, lng).ToPoint())
        loop = s2.S2Loop(inner_loops)
        loop.Normalize()  # Normalize to ensure it represents the smaller region
        polygon.InitNested([loop])  # InitNested takes a list of loops

    r_earth_km = 6371  # radius in kilometers
    print('area', polygon.GetArea())
    # Calculate the area in square kilometers
    print('area (sqkm)', polygon.GetArea() *r_earth_km ** 2)

    # Calculate the S2 cell IDs
    cell_ids = calculate_s2_cell_ids(polygon)

    print('cell_ids', cell_ids)
    # Update the s2 field (not yet)
    # obj.s2 = cell_ids
    # obj.save(update_fields=['s2'])

Editorial aside: ChatGPT did not recognize this problem, and proposed I consult maintainers. Fact is I could never have attempted this operation without its (incomplete) knowledge of the s2geometry library. And never have succeeded without its maintainers.