view pkg/models/cross.go @ 652:f5ecd1d72a6e

Cross section: Replaced naive O(N^2) segment joining with a more clever one.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Fri, 14 Sep 2018 11:44:27 +0200
parents ce18231825a2
children 7aeacd7f150b
line wrap: on
line source

package models

import (
	"bytes"
	"encoding/binary"
	"encoding/json"
	"errors"
	"fmt"
	"log"
	"math"
	"sort"
	"time"
)

type (
	GeoJSONFeature        string
	GeoJSONLineStringType string

	GeoJSONDate struct{ time.Time }

	GeoJSONLineString struct {
		Type GeoJSONLineStringType `json:"type"`
	}

	GeoJSONCoordinate struct {
		Lon float64
		Lat float64
	}

	GeoJSONCoordinateZ struct {
		Lon float64
		Lat float64
		Z   float64
	}

	CrossSectionInputProperties struct {
		Bottleneck string      `json:"bottleneck"`
		Date       GeoJSONDate `json:"date"`
	}

	GeoJSONLineCoordinates  []GeoJSONCoordinate
	GeoJSONLineCoordinatesZ []GeoJSONCoordinateZ

	GeoJSONMultiLineCoordinatesZ []GeoJSONLineCoordinatesZ

	CrossSectionInput struct {
		Type        GeoJSONFeature              `json:"type"`
		Geometry    GeoJSONLineStringType       `json:"geometry"`
		Coordinates GeoJSONLineCoordinates      `json:"coordinates"`
		Properties  CrossSectionInputProperties `json:"properties"`
	}

	CrossSectionOutput struct {
		Type        string                       `json:"type"`
		Geometry    string                       `json:"geometry"`
		Coordinates GeoJSONMultiLineCoordinatesZ `json:"coordinates"`
	}
)

var (
	errNoGeoJSONFeature        = errors.New("Not a valid GeoJSON feature")
	errNoGeoJSONLineStringType = errors.New("Not a valid GeoJSON linestring type")
	errTooLessComponents       = errors.New("Too less components in coordinate")
	errTooLessCoordinates      = errors.New("Too less coordinates")
)

func (lc *GeoJSONLineCoordinates) UnmarshalJSON(data []byte) error {
	var coords []GeoJSONCoordinate
	if err := json.Unmarshal(data, &coords); err != nil {
		return err
	}
	if len(coords) < 2 {
		return errTooLessCoordinates
	}
	*lc = GeoJSONLineCoordinates(coords)
	return nil
}

func (d *GeoJSONDate) UnmarshalJSON(data []byte) error {
	var s string
	if err := json.Unmarshal(data, &s); err != nil {
		return err
	}
	t, err := time.Parse("2006-01-02", s)
	if err != nil {
		return err
	}
	*d = GeoJSONDate{t}
	return nil
}

func (c *GeoJSONCoordinate) UnmarshalJSON(data []byte) error {
	var pos []float64
	if err := json.Unmarshal(data, &pos); err != nil {
		return err
	}
	if len(pos) < 2 {
		return errTooLessComponents
	}
	*c = GeoJSONCoordinate{Lon: pos[0], Lat: pos[1]}
	return nil
}

func (t *GeoJSONFeature) UnmarshalJSON(data []byte) error {
	var s string
	if err := json.Unmarshal(data, &s); err != nil {
		return err
	}
	if s != "Feature" {
		return errNoGeoJSONFeature
	}
	*t = GeoJSONFeature(s)
	return nil
}

func (t *GeoJSONLineStringType) UnmarshalJSON(data []byte) error {
	var s string
	if err := json.Unmarshal(data, &s); err != nil {
		return err
	}
	if s != "LineString" {
		return errNoGeoJSONLineStringType
	}
	*t = GeoJSONLineStringType(s)
	return nil
}

const (
	wkbXDR         byte   = 0
	wkbNDR         byte   = 1
	wkbPoint       uint32 = 1
	wkbPointZ      uint32 = 1000 + 1
	wkbLineString  uint32 = 2
	wkbLineStringZ uint32 = 1000 + 2
)

func (lc GeoJSONLineCoordinates) AsWKB() []byte {

	size := 1 + 4 + 4 + len(lc)*(1+4+2*8)

	buf := bytes.NewBuffer(make([]byte, 0, size))

	binary.Write(buf, binary.LittleEndian, wkbNDR)
	binary.Write(buf, binary.LittleEndian, wkbLineString)
	binary.Write(buf, binary.LittleEndian, uint32(len(lc)))

	for i := range lc {
		c := &lc[i]
		binary.Write(buf, binary.LittleEndian, math.Float64bits(c.Lon))
		binary.Write(buf, binary.LittleEndian, math.Float64bits(c.Lat))
	}

	return buf.Bytes()
}

func (cz GeoJSONCoordinateZ) MarshalJSON() ([]byte, error) {
	var buf bytes.Buffer
	fmt.Fprintf(&buf, "[%.8f,%.8f,%.8f]", cz.Lon, cz.Lat, cz.Z)
	return buf.Bytes(), nil
}

func (lcz *GeoJSONLineCoordinatesZ) Scan(src interface{}) error {
	data, ok := src.([]byte)
	if !ok {
		return errNoByteSlice
	}
	return lcz.FromWKB(data)
}

func (lcz *GeoJSONLineCoordinatesZ) FromWKB(data []byte) error {

	r := bytes.NewReader(data)

	endian, err := r.ReadByte()

	var order binary.ByteOrder

	switch {
	case err != nil:
		return err
	case endian == wkbNDR:
		order = binary.LittleEndian
	case endian == wkbXDR:
		order = binary.BigEndian
	default:
		return fmt.Errorf("unknown byte order %x", endian)
	}

	var geomType uint32
	err = binary.Read(r, order, &geomType)

	switch {
	case err != nil:
		return err
	case geomType != wkbLineStringZ:
		return fmt.Errorf("unknown geometry type %x", geomType)
	}

	var num uint32
	if err = binary.Read(r, order, &num); err != nil {
		return err
	}

	coords := make(GeoJSONLineCoordinatesZ, num)

	for i := range coords {
		c := &coords[i]
		for _, addr := range []*float64{&c.Lat, &c.Lon, &c.Z} {
			if err = binary.Read(r, order, addr); err != nil {
				return err
			}
		}
	}

	*lcz = coords

	return nil
}

func (cz GeoJSONCoordinateZ) Equals(other GeoJSONCoordinateZ) bool {
	const (
		xyEps = 1e-7
		zEps  = 1e-5
	)
	return math.Abs(cz.Lat-other.Lat) < xyEps &&
		math.Abs(cz.Lon-other.Lon) < xyEps &&
		math.Abs(cz.Z-other.Z) < zEps
}

func (cz GeoJSONCoordinateZ) mid(other GeoJSONCoordinateZ) GeoJSONCoordinateZ {
	return GeoJSONCoordinateZ{
		Lon: (cz.Lon + other.Lon) * 0.5,
		Lat: (cz.Lat + other.Lat) * 0.5,
		Z:   (cz.Z + other.Z) * 0.5,
	}
}

func deg2rad(d float64) float64 { return d * math.Pi / 180.0 }

func (cz GeoJSONCoordinateZ) Distance(other GeoJSONCoordinateZ) float64 {

	const EarthRadius = 6378137.0

	dLat := deg2rad(other.Lat - cz.Lat)
	dLng := math.Abs(deg2rad(other.Lon - cz.Lon))

	if dLng > math.Pi {
		dLng = 2*math.Pi - dLng
	}

	x := dLng * math.Cos(deg2rad((cz.Lat+other.Lat)/2.0))
	return math.Sqrt(dLat*dLat+x*x) * EarthRadius
}

func (cz GeoJSONCoordinateZ) quant() GeoJSONCoordinateZ {
	const (
		xyScale = 1e6
		zScale  = 1e5
	)
	return GeoJSONCoordinateZ{
		Lon: math.Round(cz.Lon*xyScale) / xyScale,
		Lat: math.Round(cz.Lat*xyScale) / xyScale,
		Z:   math.Round(cz.Z*zScale) / zScale,
	}
}

func (lcz GeoJSONLineCoordinatesZ) join(other GeoJSONLineCoordinatesZ) GeoJSONLineCoordinatesZ {
	l1, l2 := len(lcz), len(other)
	n := make(GeoJSONLineCoordinatesZ, l1+l2-1)
	copy(n, lcz[:l1-1])
	n[l1-1] = lcz[l1-1].mid(other[0])
	copy(n[l1:], other[1:])
	return n
}

func (ml GeoJSONMultiLineCoordinatesZ) Join() GeoJSONMultiLineCoordinatesZ {

	if len(ml) < 2 {
		return ml
	}

	start := time.Now()

	type value struct {
		coords GeoJSONLineCoordinatesZ
		order  int
	}

	heads := make(map[GeoJSONCoordinateZ]*value)

	for i, coords := range ml {
		key := coords[len(coords)-1].quant()
		if v := heads[key]; v != nil {
			delete(heads, key)
			v.coords = coords.join(v.coords)
			heads[v.coords[0].quant()] = v
		} else {
			heads[coords[0].quant()] = &value{coords, i}
		}
	}

again:
	for {
		for k, v1 := range heads {
			key := v1.coords[len(v1.coords)-1].quant()
			if k == key { // cycle detected
				continue
			}
			if v2 := heads[key]; v2 != nil {
				delete(heads, key)
				v1.coords = v1.coords.join(v2.coords)
				continue again
			}
		}
		break
	}

	// To make it deterministic sort in input order.
	tmp := make([]*value, len(heads))
	var i int
	for _, v := range heads {
		tmp[i] = v
		i++
	}
	sort.Slice(tmp, func(i, j int) bool { return tmp[i].order < tmp[j].order })

	out := make(GeoJSONMultiLineCoordinatesZ, len(tmp))
	for i, v := range tmp {
		out[i] = v.coords
	}

	log.Printf("joining took: %s\n", time.Since(start))
	log.Printf("segments before/after: %d %d\n", len(ml), len(out))

	return out
}

/*
func (ml GeoJSONMultiLineCoordinatesZ) Join() GeoJSONMultiLineCoordinatesZ {

	if len(ml) < 2 {
		return ml
	}

	start := time.Now()

	lists := make(GeoJSONMultiLineCoordinatesZ, len(ml))
	copy(lists, ml)

next:
	for j := 0; j < len(lists); {
		front := lists[j]
		tail := front[len(front)-1]

		for i := 0; i < len(lists); i++ {
			if i == j {
				continue
			}
			curr := lists[i]
			head := curr[0]
			if tail.Equals(head) {
				n := make(GeoJSONLineCoordinatesZ, len(front)+len(curr)-1)
				copy(n, front[:len(front)-1])
				n[len(front)-1] = head.mid(tail)
				copy(n[len(front):], curr[1:])
				lists[j] = n
				if i < len(lists)-1 {
					copy(lists[i:], lists[i+1:])
				}
				lists[len(lists)-1] = nil
				lists = lists[:len(lists)-1]
				continue next
			}
		}
		j++
	}

	log.Printf("joining took: %s\n", time.Since(start))
	log.Printf("segments before/after: %d %d\n", len(ml), len(lists))

	return lists
}
*/