view pkg/controllers/cross.go @ 2624:9dbaf69c7a66

Improve geoserver config to better calculate bounding boxes * Disable the use of estimated extents for the postgis storage configuration for geoserver, which is set via the gemma middleware. This way we are able to get better bounding boxes for many layers where the postgis function `ST_EstimatedExtent()` would be far off.
author Bernhard Reiter <bernhard@intevation.de>
date Wed, 13 Mar 2019 16:18:39 +0100
parents 5afca5bc1d7a
children c50d955372b9
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"
	"log"
	"net/http"
	"time"

	"gemma.intevation.de/gemma/pkg/common"
	"gemma.intevation.de/gemma/pkg/models"
	"gemma.intevation.de/gemma/pkg/octree"
)

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 octree.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(
	input interface{},
	req *http.Request,
	conn *sql.Conn,
) (jr JSONResult, err error) {

	csi := input.(*models.CrossSectionInput)

	start := time.Now()
	ctx := req.Context()

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

	log.Printf("info: loading octree took %s\n", time.Since(start))
	if err != nil {
		return
	}

	if tree == nil {
		err = 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 octree 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.Printf("info: transforming input coords took %s\n", time.Since(start))
	if err != nil {
		return
	}

	start = time.Now()

	var segments octree.MultiLineStringZ

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

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

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

		if len(line) > 0 {
			log.Printf("info: 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.Printf("info: joined length: %d\n", len(joined))
			segments = append(segments, joined...)
		}

	}
	log.Printf("info: octree traversal took %s\n", time.Since(start))

	start = time.Now()

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

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

	jr = JSONResult{
		Result: &models.CrossSectionOutput{
			Type: "Feature",
			Geometry: models.CrossSectionOutputGeometry{
				Type:        "MultiLineString",
				Coordinates: joined,
			},
			Properties: map[string]interface{}{
				"waterlevel": map[string]interface{}{
					// TODO: Fetch values from database.
					"value": float64(50),
					"when":  start.Format(common.TimeFormat),
				},
			},
		},
	}

	return
}