google / s2geometry

Computational geometry and spatial indexing on the sphere
http://s2geometry.io/
Apache License 2.0
2.29k stars 302 forks source link

trouble generating cell ids for MultiPollygons #321

Closed kgeographer closed 1 year ago

kgeographer commented 1 year ago

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)