Skip to content

trouble generating cell ids for MultiPollygons #321

@kgeographer

Description

@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 ring in geom:
            latlngs = []
            for coords in ring:
                for coord in coords:
                    lng, lat = coord
                    latlngs.append(s2.S2LatLng.FromDegrees(float(lat), float(lng)))
            loop = s2.S2Loop([latlng.ToPoint() for latlng in latlngs])

            # Create an S2Polygon from the S2Loop
            polygon = s2.S2Polygon(loop)

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

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

# 
qs = PlaceGeom.objects.filter(place__dataset='us_counties_3890')
update_s2_cell_ids(qs)

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions