Code to Convert CellUnion Back to Coordinate Polygon

[karimabedrabbo]

I had this problem where I needed to convert my list of cells back into a counterclockwise path to use in elasticsearch, so I cooked up a function to do exactly that. It definitely does not run the fastest because I only needed it for 6000 items at level 16 cells, which took roughly 3 minutes I think. Your mileage may vary. You will also need to add 1 and 2 to the Cell Union class in source. 3 is a correction to the cell union class "expand" method. 4 is the function that does it, takes as parameters a list of tokens representing your cell union, finds all the islands of the union, and creates a counter clockwise path for each island. I've only tested it for my case using standard lists of level 16 cell tokens. It should work for most other cases but take it with a grain of salt. Modify it as you need :). I went through too much pain for this, I hope this helps someone.

@classmethod
def reverse_flood_fill(cls, region, cell_union, start):
    all_nbrs = set()
    frontier = []
    all_nbrs.add(start)
    frontier.append(start)
    while len(frontier) != 0:
        cell_id = frontier.pop()
        if not region.may_intersect(Cell(cell_id)):
            continue
        yield cell_id

        neighbors = cell_id.get_edge_neighbors()
        for nbr in neighbors:
            if (not cell_union.contains(nbr)) and (nbr not in all_nbrs):
                all_nbrs.add(nbr)
                frontier.append(nbr)

@classmethod
def get_simple_reverse_covering(cls, region, cell_union, start, level):
    return cls.reverse_flood_fill(region, cell_union, CellId.from_point(start).parent(level))

    def expand(self, *args):
        if len(args) == 1 and isinstance(args[0], int):
            level = args[0]
            output = []
            level_lsb = CellId.lsb_for_level(level)
            i = self.num_cells() - 1
            while i >= 0:
                cell_id = self.__cell_ids[i]
                if cell_id.lsb() < level_lsb:
                    cell_id = cell_id.parent(level)
                    while i > 0 and cell_id.contains(self.__cell_ids[i - 1]):
                        i -= 1
                output.append(cell_id)
                output.extend(cell_id.get_all_neighbors(level))
                i -= 1
            self.__cell_ids = output
        elif len(args) == 2 and isinstance(args[0], Angle) \
                and isinstance(args[1], int):
            min_radius, max_level_diff = args
            min_level = CellId.MAX_LEVEL
            for cell_id in self.__cell_ids:
                min_level = min(min_level, cell_id.level())

            radius_level = CellId.min_width().get_max_level(min_radius.radians)
            if radius_level == 0 \
                    and min_radius.radians > CellId.min_width().get_value(0):
                self.expand(0)
            self.expand(min(min_level + max_level_diff, radius_level))
        else:
            raise NotImplementedError()

    def find_all_boundary(parent_cell_union_list, level=16):

    def find_islands(cell_union_list, level):

        tight_union = s2.CellUnion(list(map(lambda x: s2.CellId.from_token(x), cell_union_list))).denormalize(level,1)

        assert len(tight_union) > 0

        all_islands = []

        all_cells = set(map(lambda x: x.to_token(), tight_union[:]))


        while all_cells:
            all_cells_copy = all_cells.copy()

            start = all_cells_copy.pop()
            frontier = set([start])
            all_nbrs = set([start])

            while frontier:
                cell_id = frontier.pop()
                neighbors = list(map(lambda x: x.to_token(), s2.CellId.from_token(cell_id).get_edge_neighbors()))
                for nbr in neighbors:
                    if nbr in all_cells:
                        if nbr in all_cells_copy and nbr not in frontier:
                            frontier.add(nbr)
                            all_cells_copy.remove(nbr)
                            all_nbrs.add(nbr)

            all_islands.append(list(all_nbrs))

            for nbr in all_islands[-1]:
                all_cells.remove(nbr)

        return all_islands

    def findOuterBoundary(cell_union_island_list, level):

        def dist(x,y):
            return math.hypot(y[0]-x[0],y[1]-x[1])

        tight_union = s2.CellUnion(list(map(lambda x: s2.CellId.from_token(x), cell_union_island_list)))

        tight_union_expanded = s2.CellUnion(list(map(lambda x: s2.CellId.from_token(x), cell_union_island_list)))
        tight_union_expanded.expand(level)

        tight_union_bound = tight_union.get_rect_bound()

        topleft = s2.CellId.from_lat_lng(s2.LatLng.from_degrees(tight_union_bound.lat_hi().degrees, tight_union_bound.lng_lo().degrees)).parent(level)
        topleftneighborbox = s2.CellUnion(list(topleft.get_all_neighbors(level)),raw=False).get_rect_bound()


        topleft2 = s2.CellId.from_lat_lng(s2.LatLng.from_degrees(topleftneighborbox.lat_hi().degrees, topleftneighborbox.lng_lo().degrees)).parent(level)


        bottomright = s2.CellId.from_lat_lng(s2.LatLng.from_degrees(tight_union_bound.lat_lo().degrees, tight_union_bound.lng_hi().degrees)).parent(level)
        bottomrightneighborbox = s2.CellUnion(list(bottomright.get_all_neighbors(level)),raw=False).get_rect_bound()

        bottomright2 = s2.CellId.from_lat_lng(s2.LatLng.from_degrees(bottomrightneighborbox.lat_lo().degrees, bottomrightneighborbox.lng_hi().degrees)).parent(level)



        topleftpoint = s2.LatLng.from_degrees(s2.CellUnion([topleft2]).get_rect_bound().lat_hi().degrees, s2.CellUnion([topleft2]).get_rect_bound().lng_lo().degrees)
        bottomrightpoint = s2.LatLng.from_degrees(s2.CellUnion([bottomright2]).get_rect_bound().lat_lo().degrees, s2.CellUnion([bottomright2]).get_rect_bound().lng_hi().degrees)


        outerbound = s2.LatLngRect.from_point_pair(topleftpoint,bottomrightpoint)

        rc = s2.RegionCoverer()

        rc.max_level = level
        rc.min_level = level
        rc.level_mod = 1


        outercovering = s2.CellUnion(list(rc.get_simple_reverse_covering(outerbound, tight_union, topleft2.to_point(), level)))

        outercoveringexpanded = s2.CellUnion(list(rc.get_simple_reverse_covering(outerbound, tight_union, topleft2.to_point(), level)))

        outercoveringexpanded.expand(level)


        outercoveringexpanded.normalize()

        outercovering.normalize()

        tight_union_expanded.normalize()

        tight_union.normalize()


        finalintersectioninner = list(s2.CellUnion.get_intersection( outercoveringexpanded,tight_union).denormalize(level, 1))

        finalintersectionouter = list(s2.CellUnion.get_intersection( outercovering,tight_union_expanded).denormalize(level, 1))

        innertokens = set(list(map(lambda x: x.to_token(), finalintersectioninner)))
        outertokens = set(list(map(lambda x: x.to_token(), finalintersectionouter)))


        def checkccw(points):
            doubleArea = 0
            for indexp, point in enumerate(points):
                x1 = point[0]
                y1 = point[1]
                x2 = points[(indexp + 1) % len(points)][0]
                y2 = points[(indexp + 1) % len(points)][1]
                doubleArea += (x2 - x1) * (y2 + y1)
            return doubleArea

        def edgeconnections(forcoordinate_cellid_list, ccw=True, double=False):

            def algo(x):
                return (math.atan2(x[0] - mlat, x[1] - mlng) + 2 * math.pi) % (2*math.pi)

            edgeset = set()

            for thiscellid in forcoordinate_cellid_list:

                cellrep = s2.Cell(thiscellid)

                s2v0 = s2.LatLng.from_point(cellrep.get_vertex(0))
                s2v1 = s2.LatLng.from_point(cellrep.get_vertex(1))
                s2v2 = s2.LatLng.from_point(cellrep.get_vertex(2))
                s2v3 = s2.LatLng.from_point(cellrep.get_vertex(3))

                v0 = (s2v0.lat().degrees, s2v0.lng().degrees)
                v1 = (s2v1.lat().degrees, s2v1.lng().degrees)
                v2 = (s2v2.lat().degrees, s2v2.lng().degrees)
                v3 = (s2v3.lat().degrees, s2v3.lng().degrees)

                s2center = s2.LatLng.from_point(cellrep.get_center())

                mlat = s2center.lat().degrees
                mlng = s2center.lng().degrees

                edgeset.add((v0, v1))
                edgeset.add((v1, v0))

                edgeset.add((v0, v3))
                edgeset.add((v3, v0))

                edgeset.add((v1, v2))
                edgeset.add((v2, v1))

                edgeset.add((v3, v2))
                edgeset.add((v2, v3))

            return edgeset


        a = edgeconnections(finalintersectioninner,double=True,ccw=False)
        b = edgeconnections(finalintersectionouter,double=True,ccw=True)


        boundingedgeset = a.intersection(b)




        specificedge = boundingedgeset.pop()
        finalpath = [specificedge[0]]
        seencoords = set([specificedge[0]])
        boundingedgeset.add(specificedge)
        fresh = True

        edgecoorddict = {}
        for ttedge in boundingedgeset:
            if (ttedge[0][0], ttedge[0][1]) in edgecoorddict:
                edgecoorddict[(ttedge[0][0], ttedge[0][1])].add(ttedge)
            else:
                edgecoorddict[(ttedge[0][0], ttedge[0][1])] = set([ttedge])
            if (ttedge[1][0], ttedge[1][1]) in edgecoorddict:
                edgecoorddict[(ttedge[1][0], ttedge[1][1])].add(ttedge)
            else:
                edgecoorddict[(ttedge[1][0], ttedge[1][1])] = set([ttedge])


        while finalpath[0] != specificedge[0] or fresh:
            fresh = False

            boundingedgeset.remove(specificedge)

            #start of edge
            if specificedge[0] not in seencoords:
                finalpath.append(specificedge[0])
                seencoords.add(specificedge[0])

            coorspcoordset = edgecoorddict[(specificedge[1][0],specificedge[1][1])]

            copyset = coorspcoordset.copy()

            copyset.discard(specificedge)

            copyset.discard((specificedge[1], specificedge[0]))

            if len(copyset) == 1:
                specificedge = copyset.pop()
            elif len(copyset) == 2:
                cand1, cand2 = copyset.pop(), copyset.pop()
                if not cand1[1] == specificedge[1]:
                    specificedge = cand1
                elif not cand2[1] == specificedge[1]:
                    specificedge = cand2


        actualfinalpath = finalpath + [finalpath[0]]
        ccwarea = checkccw(actualfinalpath)


        if ccwarea <= 0:
            return finalintersectioninner, finalintersectionouter, actualfinalpath
        else:
            return finalintersectioninner, finalintersectionouter, actualfinalpath[::-1]

    def reformatToSingleBoundaryOSM(boundary):
        thestr = ''
        thestr += '(('
        for thetuple in boundary:
            thestr += str(thetuple[1]) + ' ' + str(thetuple[0]) + ','
        thestr = thestr[:-1]
        thestr += ')),'
        thestr = thestr[:-1]
        return "MULTIPOLYGON(" + thestr + ")"

    islandreturnlist = []

    for island in find_islands(parent_cell_union_list, level):
        cellboundaryinner, cellboundaryouter, ccwpath = findOuterBoundary(island, level)

        osmccwpath = reformatToSingleBoundaryOSM(ccwpath)


        newdict = {}

        newdict["island"] = island
        newdict["cellboundaryinner"] = list(map(lambda x: x.to_token(), cellboundaryinner))
        newdict["cellboundaryouter"] = list(map(lambda x: x.to_token(), cellboundaryouter))
        newdict["ccwpath"] = ccwpath
        newdict["osmccwpath"] = osmccwpath

        islandreturnlist.append(newdict)

    return islandreturnlist