view pkg/imports/wkb.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 41a67619c170
children 6270951dda28
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 imports

import (
	"bytes"
	"encoding/binary"
	"fmt"
	"math"

	shp "github.com/jonas-p/go-shp"

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

type (
	pointSlice        []float64
	lineSlice         [][]float64
	multiLineSlice    []lineSlice
	polygonSlice      [][][]float64
	multiPolygonSlice []polygonSlice
)

func newPointFeature(newProperties func() interface{}) func() (string, interface{}) {
	return func() (string, interface{}) { return "Point", newProperties() }
}

func newMultiLineFeature(
	newProperties func() interface{},
) func() (string, interface{}) {
	return func() (string, interface{}) {
		return "MultiLineString", newProperties()
	}
}

func (ls lineSlice) toWKB(buf *bytes.Buffer) {
	binary.Write(buf, binary.LittleEndian, wkb.NDR)
	binary.Write(buf, binary.LittleEndian, wkb.LineString)
	binary.Write(buf, binary.LittleEndian, uint32(len(ls)))

	for _, c := range ls {
		var lat, lon float64
		if len(c) > 0 {
			lat = c[0]
		}
		if len(c) > 1 {
			lon = c[1]
		}
		binary.Write(buf, binary.LittleEndian, math.Float64bits(lat))
		binary.Write(buf, binary.LittleEndian, math.Float64bits(lon))
	}
}

func (ls lineSlice) asWKB() []byte {

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

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

	return buf.Bytes()
}

func (ls lineSlice) LinearRingGeom() wkb.LinearRingGeom {
	lr := make(wkb.LinearRingGeom, len(ls))
	for i, v := range ls {
		lr[i].X = v[0]
		lr[i].Y = v[1]
	}
	return lr
}

func (mls multiLineSlice) asWKB() []byte {

	size := 1 + 4 + 4

	for _, ls := range mls {
		size += 1 + 4 + 4 + len(ls)*(2*8)
	}

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

	binary.Write(buf, binary.LittleEndian, wkb.NDR)
	binary.Write(buf, binary.LittleEndian, wkb.MultiLineString)
	binary.Write(buf, binary.LittleEndian, uint32(len(mls)))

	for _, ls := range mls {
		ls.toWKB(buf)
	}

	return buf.Bytes()
}

func (p pointSlice) asWKB() []byte {

	size := 1 + 4 + 2*8

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

	binary.Write(buf, binary.LittleEndian, wkb.NDR)
	binary.Write(buf, binary.LittleEndian, wkb.Point)

	var lat, lon float64
	if len(p) > 0 {
		lat = p[0]
	}
	if len(p) > 1 {
		lon = p[1]
	}
	binary.Write(buf, binary.LittleEndian, math.Float64bits(lat))
	binary.Write(buf, binary.LittleEndian, math.Float64bits(lon))

	return buf.Bytes()
}

func (ps polygonSlice) asWKB() []byte {
	if ps == nil {
		return nil
	}
	// pre-calculate size to avoid reallocations.
	size := 1 + 4 + 4
	for _, ring := range ps {
		size += 4 + len(ring)*2*8
	}

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

	binary.Write(buf, binary.LittleEndian, wkb.NDR)
	binary.Write(buf, binary.LittleEndian, wkb.Polygon)
	binary.Write(buf, binary.LittleEndian, uint32(len(ps)))

	for _, ring := range ps {
		binary.Write(buf, binary.LittleEndian, uint32(len(ring)))
		for _, v := range ring {
			var lat, lon float64
			if len(v) > 0 {
				lat = v[0]
			}
			if len(v) > 1 {
				lon = v[1]
			}
			binary.Write(buf, binary.LittleEndian, math.Float64bits(lat))
			binary.Write(buf, binary.LittleEndian, math.Float64bits(lon))
		}
	}

	return buf.Bytes()
}

func shapeToPolygon(s shp.Shape) (polygonSlice, error) {
	switch p := s.(type) {
	case *shp.Polygon:
		return toPolygon(p.NumParts, p.Parts, p.Points), nil
	case *shp.PolygonZ:
		return toPolygon(p.NumParts, p.Parts, p.Points), nil
	case *shp.PolygonM:
		return toPolygon(p.NumParts, p.Parts, p.Points), nil
	}
	return nil, fmt.Errorf("unsupported shape type %T", s)
}

func toPolygon(numParts int32, parts []int32, points []shp.Point) polygonSlice {
	out := make(polygonSlice, numParts)
	var pos int32

	for i := range out {
		var howMany int32
		if i+1 >= len(parts) {
			howMany = int32(len(points)) - pos
		} else {
			howMany = parts[i+1] - parts[i]
		}

		line := make([][]float64, howMany)
		vertices := make([]float64, 2*howMany)
		for j := int32(0); j < howMany; j, pos = j+1, pos+1 {
			p := &points[pos]
			vertex := vertices[j*2 : j*2+2]
			vertex[0], vertex[1] = p.X, p.Y
			line[j] = vertex
		}
		out[i] = line
	}
	return out
}

func (ps polygonSlice) MultiPolygonGeom() wkb.MultiPolygonGeom {

	var mp wkb.MultiPolygonGeom
	var curr wkb.PolygonGeom

	for _, r := range ps {
		lr := lineSlice(r).LinearRingGeom()
		// A counter clockwise ring starts a new polygon.
		if lr.CCW() {
			if len(curr) > 0 {
				mp = append(mp, curr)
				curr = wkb.PolygonGeom{}
			}
		}
		curr = append(curr, lr)
	}

	if len(curr) > 0 {
		mp = append(mp, curr)
	}

	return mp
}