view pkg/models/cross.go @ 5520:05db984d3db1

Improve performance of bottleneck area calculation Avoid buffer calculations by replacing them with simple distance comparisons and calculate the boundary of the result geometry only once per iteration. In some edge cases with very large numbers of iterations, this reduced the runtime of a bottleneck import by a factor of more than twenty.
author Tom Gottfried <tom@intevation.de>
date Thu, 21 Oct 2021 19:50:39 +0200
parents 4847ac70103a
children 1222b777f51f
line wrap: on
line source

// This is Free Software under GNU Affero General Public License v >= 3.0
// without warranty, see README.md and license for details.
//
// SPDX-License-Identifier: AGPL-3.0-or-later
// License-Filename: LICENSES/AGPL-3.0.txt
//
// Copyright (C) 2018 by via donau
//   – Österreichische Wasserstraßen-Gesellschaft mbH
// Software engineering by Intevation GmbH
//
// Author(s):
//  * Sascha L. Teichmann <sascha.teichmann@intevation.de>

package models

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

	"gemma.intevation.de/gemma/pkg/common"
	"gemma.intevation.de/gemma/pkg/wkb"
)

type (
	GeoJSONFeature        string
	GeoJSONLineStringType string

	GeoJSONDate struct{ time.Time }

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

	GeoJSONCoordinate struct {
		Lat float64
		Lon float64
	}

	GeoJSONCoordinateZ struct {
		Lat float64
		Lon float64
		Z   float64
	}

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

	GeoJSONLineCoordinates  []GeoJSONCoordinate
	GeoJSONLineCoordinatesZ []GeoJSONCoordinateZ

	GeoJSONMultiLineCoordinatesZ []GeoJSONLineCoordinatesZ

	CrossSectionInputGeometry struct {
		Type        GeoJSONLineStringType  `json:"type"`
		Coordinates GeoJSONLineCoordinates `json:"coordinates"`
	}

	CrossSectionInput struct {
		Type       GeoJSONFeature              `json:"type"`
		Geometry   CrossSectionInputGeometry   `json:"geometry"`
		Properties CrossSectionInputProperties `json:"properties"`
	}

	CrossSectionOutputGeometry struct {
		Type        string                       `json:"type"`
		Coordinates GeoJSONMultiLineCoordinatesZ `json:"coordinates"`
	}

	CrossSectionOutput struct {
		Type       string                     `json:"type"`
		Geometry   CrossSectionOutputGeometry `json:"geometry"`
		Properties map[string]interface{}     `json:"properties"`
	}
)

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(common.DateFormat, 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{Lat: pos[0], Lon: 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
}

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

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

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

	binary.Write(buf, binary.LittleEndian, wkb.NDR)
	binary.Write(buf, binary.LittleEndian, wkb.LineString)
	binary.Write(buf, binary.LittleEndian, uint32(len(lc)))

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

	return buf.Bytes()
}

func (cz GeoJSONCoordinateZ) MarshalJSON() ([]byte, error) {
	var buf bytes.Buffer
	fmt.Fprintf(&buf, "[%.8f,%.8f,%.8f]", cz.Lat, cz.Lon, 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 == wkb.NDR:
		order = binary.LittleEndian
	case endian == wkb.XDR:
		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 != wkb.LineStringZ:
		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 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 (mls *GeoJSONMultiLineCoordinatesZ) FromWKB(data []byte) error {
	r := bytes.NewReader(data)

	var order binary.ByteOrder

	endian, err := r.ReadByte()

	switch {
	case err != nil:
		return err
	case endian == wkb.NDR:
		order = binary.LittleEndian
	case endian == wkb.XDR:
		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 != wkb.MultiLineStringZ:
		return fmt.Errorf("unknown geometry type %d", geomType)
	}

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

	lines := make(GeoJSONMultiLineCoordinatesZ, numLines)

	for i := range lines {
		endian, err = r.ReadByte()

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

		err = binary.Read(r, order, &geomType)
		switch {
		case err != nil:
			return err
		case geomType != wkb.LineStringZ:
			return fmt.Errorf("unknown geometry type %d", geomType)
		}

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

		points := make(GeoJSONLineCoordinatesZ, numPoints)
		for j := range points {
			var lat, lon, z uint64
			if err = binary.Read(r, order, &lat); err != nil {
				return err
			}
			if err = binary.Read(r, order, &lon); err != nil {
				return err
			}
			if err = binary.Read(r, order, &z); err != nil {
				return err
			}
			c := &points[j]
			c.Lat = math.Float64frombits(lat)
			c.Lon = math.Float64frombits(lon)
			c.Z = math.Float64frombits(z)
		}
		lines[i] = points
	}

	*mls = lines

	return nil
}

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