view pkg/controllers/diff.go @ 4574:daed8a92024a iso-areas

Generate iso ares for sounding differences, too.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Mon, 07 Oct 2019 13:13:14 +0200
parents 4ca884dfc470
children c657dec6b0fa
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, 2019 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 (
	"database/sql"
	"fmt"
	"log"
	"net/http"
	"sync"
	"time"

	"golang.org/x/sync/semaphore"

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

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

const (
	contourTolerance = 0.1
	contourStep      = 0.1
)

const (
	// isoCellSize is the side length of a raster cell when tracing
	// iso areas.
	isoCellSize = 0.5
)

const (
	diffIDSQL = `
SELECT sd.id FROM
  caching.sounding_differences sd JOIN
  waterway.sounding_results srm ON sd.minuend = srm.id JOIN
  waterway.sounding_results srs ON sd.subtrahend = srs.id
WHERE srm.bottleneck_id = srs.bottleneck_id AND
  srm.bottleneck_id = $1 AND
  srm.date_info = $2::date AND
  srs.date_info = $3::date
`
	insertDiffSQL = `
WITH soundings AS (
  SELECT sr.id AS id, sr.date_info AS date_info FROM
  waterway.sounding_results sr
  WHERE sr.bottleneck_id = $1
)
INSERT INTO caching.sounding_differences (minuend, subtrahend)
SELECT m.id, s.id FROM soundings m, soundings s
WHERE m.date_info = $2::date AND s.date_info = $3::date
RETURNING id
`
	insertDiffContourSQL = `
INSERT INTO caching.sounding_differences_contour_lines (
  sounding_differences_id,
  height,
  lines
)
SELECT
  $5,
  $4,
  ST_Transform(
    ST_Multi(
      ST_CollectionExtract(
        ST_SimplifyPreserveTopology(
          ST_Multi(ST_Collectionextract(
            ST_MakeValid(ST_GeomFromWKB($1, $2::integer)), 2)),
          $3
        ),
        2
      )
    ),
    4326
  )
`
	insertDiffIsoAreasQL = `
INSERT INTO caching.sounding_differences_iso_areas (
  sounding_differences_id,
  height,
  areas
)
SELECT
  $1,
  $2,
  ST_Transform(
    ST_Multi(
      ST_CollectionExtract(
        ST_SimplifyPreserveTopology(
          ST_Multi(ST_Collectionextract(
            ST_MakeValid(ST_GeomFromWKB($4, $3::integer)), 3)),
          $5
        ),
        3
      )
    ),
    4326
  )
`
)

// Only allow three diffence calculation at once.
// TODO: Make this configurable?
var diffCalculationSemaphore = semaphore.NewWeighted(int64(3))

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

	begin := time.Now()
	start := begin

	dci := mw.JSONInput(req).(*models.DiffCalculationInput)

	ctx := req.Context()

	conn := mw.JSONConn(req)

	var id int64
	err = conn.QueryRowContext(
		ctx,
		diffIDSQL,
		dci.Bottleneck,
		dci.Minuend.Time,
		dci.Subtrahend.Time,
	).Scan(&id)

	switch {
	case err == sql.ErrNoRows:
		// We need to calculate it.
	case err != nil:
		return
	default:
		// We already have this diff
		jr = mw.JSONResult{
			Result: map[string]int64{"id": id},
		}
		return
	}

	// DoS counter measure.
	if err = diffCalculationSemaphore.Acquire(ctx, 1); err != nil {
		return
	}
	defer diffCalculationSemaphore.Release(1)

	minuendTree, err := octree.FromCache(
		ctx, conn,
		dci.Bottleneck, dci.Minuend.Time)

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

	if minuendTree == nil {
		err = mw.JSONError{
			Code: http.StatusNotFound,
			Message: fmt.Sprintf("Cannot find survey for %s/%s.",
				dci.Bottleneck,
				dci.Minuend.Format(common.DateFormat)),
		}
		return
	}

	start = time.Now()

	subtrahendTree, err := octree.FromCache(
		ctx, conn,
		dci.Bottleneck, dci.Subtrahend.Time)

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

	if subtrahendTree == nil {
		err = mw.JSONError{
			Code: http.StatusNotFound,
			Message: fmt.Sprintf("Cannot find survey for %s/%s.",
				dci.Bottleneck,
				dci.Subtrahend.Format(common.DateFormat)),
		}
		return
	}

	// We need a slow path implementation for this.
	if minuendTree.EPSG != subtrahendTree.EPSG {
		err = mw.JSONError{
			Code: http.StatusInternalServerError,
			Message: "Calculating differences between two different " +
				"EPSG code octrees are not supported, yet.",
		}
		return
	}

	start = time.Now()
	points := minuendTree.Diff(subtrahendTree)
	log.Printf("info: A - B took %v\n", time.Since(start))

	// The Triangulation and the loading of the clipping
	// polygon can be done concurrently.

	jobs := make(chan func())

	wg := new(sync.WaitGroup)
	for i := 0; i < 2; i++ {
		wg.Add(1)
		go func() {
			defer wg.Done()
			for job := range jobs {
				job()
			}
		}()
	}

	var (
		tri     *octree.Triangulation
		triErr  error
		clip    *octree.Polygon
		clipErr error
	)

	jobs <- func() {
		start := time.Now()
		tri, triErr = points.Triangulate()
		log.Printf("info: triangulation took %v\n", time.Since(start))
	}

	jobs <- func() {
		start := time.Now()
		clip, clipErr = octree.LoadClippingPolygon(
			ctx, conn,
			minuendTree.EPSG,
			dci.Bottleneck,
			dci.Minuend.Time,
			dci.Subtrahend.Time)
		log.Printf("info: loading clipping polygon took %v\n", time.Since(start))
	}
	close(jobs)
	wg.Wait()

	switch {
	case triErr != nil && clipErr != nil:
		err = fmt.Errorf("%v %v", triErr, clipErr)
		return
	case triErr != nil:
		err = triErr
		return
	case clipErr != nil:
		err = clipErr
		return
	}

	start = time.Now()
	tin := tri.Tin()
	removed := tin.Clip(clip)
	log.Printf("info: clipping TIN took %v\n", time.Since(start))

	start = time.Now()
	log.Printf("info: Number of triangles to clip: %d\n", len(removed))
	builder := octree.NewBuilder(tin)
	builder.Build(removed)
	log.Printf("info: building octree took %v\n", time.Since(start))

	tree := builder.Tree()

	log.Printf("info: min/max: %f %f\n", tree.Min.Z, tree.Max.Z)

	start = time.Now()

	// XXX: Maybe we should start this transaction earlier!?
	var tx *sql.Tx
	if tx, err = conn.BeginTx(ctx, nil); err != nil {
		return
	}
	defer tx.Rollback()

	var heights []float64

	heights, err = octree.LoadClassBreaks(
		ctx, tx,
		"morphology_classbreaks_compare")
	if err != nil {
		log.Printf("warn: Loading class breaks failed: %v\n", err)
		heights = octree.SampleDiffHeights(
			tree.Min.Z, tree.Max.Z, contourStep)
	} else {
		heights = octree.ExtrapolateClassBreaks(heights, tree.Min.Z, tree.Max.Z)
		// heights = octree.InBetweenClassBreaks(heights, 0.05, 2)
	}

	log.Printf("info: num heights: %d\n", len(heights))

	var stmt *sql.Stmt
	if stmt, err = tx.PrepareContext(ctx, insertDiffContourSQL); err != nil {
		return
	}
	defer stmt.Close()

	var isoStmt *sql.Stmt
	if isoStmt, err = tx.PrepareContext(ctx, insertDiffIsoAreasQL); err != nil {
		return
	}
	defer isoStmt.Close()

	if err = tx.QueryRowContext(
		ctx,
		insertDiffSQL,
		dci.Bottleneck,
		dci.Minuend.Time,
		dci.Subtrahend.Time,
	).Scan(&id); err != nil {
		return
	}

	heights = common.DedupFloat64s(heights)

	octree.DoContours(tree, heights, func(res *octree.ContourResult) {
		if err == nil && len(res.Lines) > 0 {
			_, err = stmt.ExecContext(
				ctx,
				res.Lines.AsWKB2D(),
				minuendTree.EPSG,
				contourTolerance,
				res.Height,
				id,
			)
		}
	})

	areas := tree.TraceAreas(heights, isoCellSize)

	for i, a := range areas {
		if len(a) == 0 {
			continue
		}
		if _, err = isoStmt.ExecContext(
			ctx,
			id, heights[i], minuendTree.EPSG,
			a.AsWKB(),
			contourTolerance,
		); err != nil {
			return
		}
	}

	log.Printf("info: calculating and storing iso lines took %v\n",
		time.Since(start))

	if err != nil {
		return
	}

	if err = tx.Commit(); err != nil {
		log.Printf("info: difference calculation failed after %v\n",
			time.Since(begin))
		return
	}
	log.Printf("info: difference calculation succeed after %v\n",
		time.Since(begin))

	jr = mw.JSONResult{
		Result: map[string]int64{"id": id},
	}
	return
}