I'm trying to get the largest possible cells that cover but do not exceed a multi-polygon. I've tried a few different approaches to solve this but no matter what I do I end up with the same CellIDs for each loop (polygon).
package main
import (
"encoding/json"
"fmt"
"os"
"github.com/golang/geo/s2"
"github.com/twpayne/go-geom"
"github.com/twpayne/go-geom/encoding/geojson"
)
type Location struct {
Geometry json.RawMessage `json:"geometry"`
}
func main() {
f, err := os.Open("./geo.json")
if err != nil {
panic(err)
}
s, _ := f.Stat()
b := make([]byte, s.Size())
f.Read(b)
loc := &Location{}
var g geom.T
err = json.Unmarshal(b, &loc)
if err != nil {
panic(err)
}
err = geojson.Unmarshal(loc.Geometry, &g)
if err != nil {
panic(err)
}
mp := g.(*geom.MultiPolygon)
coords := mp.Coords()
for _, polygon := range coords {
var xv float64
var xi = -1
points := []s2.Point{}
for _, l := range polygon {
fmt.Println(len(l))
for i := 0; i < len(l)-1; i++ {
p := l[i]
x := p.X()
latlng := s2.LatLngFromDegrees(p.Y(), x)
point := s2.PointFromLatLng(latlng)
points = append(points, point)
if xi == -1 || x > xv {
xv = x
xi = i
}
}
}
if xi == -1 {
continue
}
c := s2.NewRegionCoverer()
// changing these values has no impact on the output
c.MaxLevel = 30
c.MinLevel = 1
points = append(points[xi:], points[:xi]...)
loop := s2.LoopFromPoints(points)
if loop.Validate() != nil {
fmt.Println(loop.Validate())
}
if loop.IsHole() {
panic("loop is hole")
}
u := c.Covering(loop)
printUnionBound(u.CellUnionBound())
// this produces the same output as using the RegionCoverer
// printUnionBound(loop.CellUnionBound())
// polygon := s2.PolygonFromLoops([]*s2.Loop{loop})
// printUnionBound(polygon.CellUnionBound())
}
}
func printUnionBound(cellIDs []s2.CellID) {
l := len(cellIDs)
for i := 0; i < l; i++ {
cellID := cellIDs[i]
if i == l-1 {
fmt.Printf("%d\n", cellID)
} else {
fmt.Printf("%d, ", cellID)
}
}
}
I'm trying to get the largest possible cells that cover but do not exceed a multi-polygon. I've tried a few different approaches to solve this but no matter what I do I end up with the same CellIDs for each loop (polygon).
The code is:
I'm not sure what I'm doing wrong.
The geojson that I'm experimenting with is:
which, when loaded into QGIS produces:
I'm not accounting for holes at the moment as this geojson doesn't have any.