view pkg/controllers/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 5f47eeea988d
children
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 controllers

import (
	"context"
	"database/sql"
	"fmt"
	"net/http"
	"time"

	"gemma.intevation.de/gemma/pkg/log"
	"gemma.intevation.de/gemma/pkg/mesh"
	"gemma.intevation.de/gemma/pkg/models"

	mw "gemma.intevation.de/gemma/pkg/middleware"
)

func reproject(
	ctx context.Context,
	rp *models.Reprojector,
	src models.GeoJSONLineCoordinates,
) (models.GeoJSONLineCoordinates, error) {

	dst := make(models.GeoJSONLineCoordinates, len(src))
	for i, s := range src {
		var err error
		if dst[i].Lat, dst[i].Lon, err = rp.Reproject(
			ctx,
			s.Lat, s.Lon,
		); err != nil {
			return nil, err
		}
	}
	return dst, nil
}

const projectBackSQL = `
SELECT ST_AsBinary(
  ST_Transform(ST_GeomFromWKB($2, $1::integer), 4326))`

func projectBack(
	ctx context.Context,
	line mesh.MultiLineStringZ,
	epsg uint32,
	conn *sql.Conn,
) (models.GeoJSONMultiLineCoordinatesZ, error) {

	var mls models.GeoJSONMultiLineCoordinatesZ
	err := conn.QueryRowContext(
		ctx, projectBackSQL,
		epsg, line.AsWKB(),
	).Scan(&mls)

	return mls, err
}

func crossSection(req *http.Request) (jr mw.JSONResult, err error) {

	csi := mw.JSONInput(req).(*models.CrossSectionInput)

	start := time.Now()
	ctx := req.Context()
	conn := mw.JSONConn(req)

	tree, err := mesh.FromCache(
		ctx, conn,
		csi.Properties.Bottleneck, csi.Properties.Date.Time)

	log.Infof("loading mesh took %s\n", time.Since(start))
	if err != nil {
		return
	}

	if tree == nil {
		err = mw.JSONError{
			Code: http.StatusNotFound,
			Message: fmt.Sprintf("Cannot find survey for %s/%s.",
				csi.Properties.Bottleneck,
				csi.Properties.Date.Time),
		}
		return
	}

	// The coordinate system of the mesh is an UTM projection.
	// The input coordinates are in WGS84.
	// So we need to reproject them.

	start = time.Now()

	var rp *models.Reprojector
	if rp, err = models.NewReprojector(
		ctx, conn,
		models.WGS84, tree.EPSG(),
	); err != nil {
		return
	}
	defer rp.Close()

	coords, err := reproject(ctx, rp, csi.Geometry.Coordinates)

	log.Infof("transforming input coords took %s\n", time.Since(start))
	if err != nil {
		return
	}

	start = time.Now()

	var segments mesh.MultiLineStringZ

	for i := 0; i < len(coords)-1; i++ {
		c1 := &coords[i]
		c2 := &coords[i+1]

		verticalLine := mesh.NewVerticalLine(c1.Lat, c1.Lon, c2.Lat, c2.Lon)

		var line mesh.MultiLineStringZ
		tree.Vertical(c1.Lat, c1.Lon, c2.Lat, c2.Lon, func(t *mesh.Triangle) {
			if ls := verticalLine.Intersection(t); len(ls) > 0 {
				line = append(line, ls)
			}
		})

		if len(line) > 0 {
			log.Infof("line length: %d\n", len(line))
			// They are all on the segment (c1.Lat, c1.Lon) - (c2.Lat, c2.Lon).
			// Sort them by project them on this line.
			joined := line.JoinOnLine(c1.Lat, c1.Lon, c2.Lat, c2.Lon)
			log.Infof("joined length: %d\n", len(joined))
			segments = append(segments, joined...)
		}

	}
	log.Infof("mesh traversal took %s\n", time.Since(start))

	start = time.Now()

	var joined models.GeoJSONMultiLineCoordinatesZ
	joined, err = projectBack(
		ctx,
		segments, tree.EPSG(),
		conn,
	)

	log.Infof("projecting back took %s\n", time.Since(start))
	if err != nil {
		return
	}

	jr = mw.JSONResult{
		Result: &models.CrossSectionOutput{
			Type: "Feature",
			Geometry: models.CrossSectionOutputGeometry{
				Type:        "MultiLineString",
				Coordinates: joined,
			},
		},
	}

	return
}