Convert golang S2 geometries to Simple Features

conversions.go

package spatial

import (
	"fmt"
	"strconv"

	"github.com/golang/geo/s2"
	"github.com/peterstace/simplefeatures/geom"
)

const (
	ORIENTATION_UNKNOWN uint8 = iota
	ORIENTATION_COUNTER_CLOCKWISE
	ORIENTATION_CLOCKWISE
)

func SFPolygonToS2Polygon(geometry geom.Geometry) *s2.Polygon {
	polygon := geometry.MustAsPolygon()
	num := polygon.NumInteriorRings()
	loops := make([]*s2.Loop, 1+num)

	shell := polygon.ExteriorRing()
	loops[0] = SFLineStringToS2Loop(&shell)

	for i := 0; i < num; i++ {
		hole := polygon.InteriorRingN(i)
		loops[1+i] = SFLineStringToS2Loop(&hole)
	}

	return s2.PolygonFromOrientedLoops(loops)
}

// https://s2geometry.io/devguide/basic_types#s2loop
func SFLineStringToS2Loop(linestring *geom.LineString) *s2.Loop {
	orient := orientation(linestring)

	seq := linestring.Coordinates()
	len := seq.Length()

	// s2: "It consists of a single chain of vertices where the first
	// vertex is implicitly connected to the last.""
	if linestring.IsClosed() {
		len -= 1
	}

	points := make([]s2.Point, len)

	// s2: "All loops are defined to have a CCW orientation"
	for s := 0; s < len; s++ {
		if orient == ORIENTATION_COUNTER_CLOCKWISE {
			points[s] = SFCoordinatesToS2Point(seq.Get(s))
		} else {
			points[len-1-s] = SFCoordinatesToS2Point(seq.Get(s))
		}
	}

	return s2.LoopFromPoints(points)
}

func SFLineStringToS2Polyline(geometry geom.Geometry) *s2.Polyline {
	linestring := geometry.MustAsLineString()
	seq := linestring.Coordinates()
	vertices := make([]s2.LatLng, seq.Length())

	for s := 0; s < seq.Length(); s++ {
		vertices[s] = SFCoordinatesToS2LatLng(seq.Get(s))
	}

	return s2.PolylineFromLatLngs(vertices)
}

func SFCoordinatesToS2Point(coords geom.Coordinates) s2.Point {
	return s2.PointFromLatLng(SFCoordinatesToS2LatLng(coords))
}
func SFCoordinatesToS2LatLng(coords geom.Coordinates) s2.LatLng {
	return s2.LatLngFromDegrees(coords.Y, coords.X)
}

func SFPointToS2Cell(point geom.Point) s2.Cell {
	coords, _ := point.Coordinates()
	return s2.CellFromLatLng(SFCoordinatesToS2LatLng(coords))
}

func S2CellIDToSFGeometry(cellID s2.CellID) geom.Geometry {
	return S2CellToSFGeometry(s2.CellFromCellID(cellID))
}

// s2: Vertex returns the k-th vertex of the cell (k = 0,1,2,3) in CCW order
// s2: (lower left, lower right, upper right, upper left in the UV plane).
// @todo: any advantage to converting to a polygon?
func S2CellToSFGeometry(cell s2.Cell) geom.Geometry {
	v0 := s2.LatLngFromPoint(cell.Vertex(0))
	v1 := s2.LatLngFromPoint(cell.Vertex(1))
	v2 := s2.LatLngFromPoint(cell.Vertex(2))
	v3 := s2.LatLngFromPoint(cell.Vertex(3))

	floats := []float64{
		v0.Lng.Degrees(),
		v0.Lat.Degrees(),
		v1.Lng.Degrees(),
		v1.Lat.Degrees(),
		v2.Lng.Degrees(),
		v2.Lat.Degrees(),
		v3.Lng.Degrees(),
		v3.Lat.Degrees(),
		v0.Lng.Degrees(),
		v0.Lat.Degrees(),
	}

	shell, _ := geom.NewLineString(
		geom.NewSequence(floats, geom.DimXY),
		geom.DisableAllValidations,
	)

	polygon, _ := geom.NewPolygon(
		[]geom.LineString{shell},
		geom.DisableAllValidations,
	)

	return polygon.AsGeometry()
}

// @todo: can this be simplified?
func SFPointFromLatLng(lat float64, lon float64) geom.Geometry {
	point, _ := geom.NewPoint(
		geom.Coordinates{XY: geom.XY{X: lon, Y: lat}},
		geom.DisableAllValidations,
	)
	return point.AsGeometry()
}

func SFPointFromLatLngStrings(latitude string, longitude string) (geom.Geometry, error) {
	lat, err := strconv.ParseFloat(latitude, 64)
	if err != nil || (lat < -90.0 || lat > 90.0) {
		return geom.Geometry{}, fmt.Errorf("invalid latitude")
	}

	lon, err := strconv.ParseFloat(longitude, 64)
	if err != nil || (lon < -180.0 || lon > 180.0) {
		return geom.Geometry{}, fmt.Errorf("invalid longitude")
	}

	return SFPointFromLatLng(lat, lon), nil
}

// Discover winding order
// https://github.com/peterstace/simplefeatures/blob/560fd6d9e24316df1ee05c2fab00916a179f8a85/geom/type_polygon.go#L391
func orientation(lr *geom.LineString) uint8 {
	// This is the "Shoelace Formula".
	var sum float64
	seq := lr.Coordinates()
	n := seq.Length()
	for i := 0; i < n; i++ {
		pt0 := seq.GetXY(i)
		pt1 := seq.GetXY((i + 1) % n)
		sum += (pt1.X + pt0.X) * (pt1.Y - pt0.Y)
	}
	ret := sum / 2

	if ret < 0 {
		return ORIENTATION_CLOCKWISE
	} else if ret > 0 {
		return ORIENTATION_COUNTER_CLOCKWISE
	}

	return ORIENTATION_UNKNOWN
}