view pkg/imports/sr.go @ 1317:5443f5c9154c

Added missing authors names in Go files.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Sat, 24 Nov 2018 10:32:06 +0100
parents 5aeda02c51b9
children d753ce6cf588
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>
//  * Tom Gottfried <tom.gottfried@intevation.de>

package imports

import (
	"archive/zip"
	"bufio"
	"context"
	"crypto/sha1"
	"database/sql"
	"encoding/hex"
	"encoding/json"
	"errors"
	"fmt"
	"io"
	"math"
	"os"
	"path"
	"path/filepath"
	"strconv"
	"strings"
	"time"

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

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

type SoundingResult struct {
	Dir string `json:"dir"`

	// Override data
	Date           *models.SoundingResultDate `json:"date,omitempty"`
	Bottleneck     *string                    `json:"bottleneck,omitempty"`
	EPSG           *uint                      `json:"epsg,omitempty"`
	DepthReference *string                    `json:"depth-reference,omitempty"`
}

const (
	contourStepWidth = 0.1
	contourTolerance = 0.1
)

const SRJobKind JobKind = "sr"

type srJobCreator struct{}

func init() {
	RegisterJobCreator(SRJobKind, srJobCreator{})
}

func (srJobCreator) Create(_ JobKind, data string) (Job, error) {
	sr := new(SoundingResult)
	if err := sr.FromString(data); err != nil {
		return nil, err
	}
	return sr, nil
}

func (srJobCreator) Depends() []string {
	return []string{
		"waterway.sounding_results",
		"waterway.sounding_results_contour_lines",
		"waterway.bottlenecks",
	}
}

func (srJobCreator) StageDone(
	tx *sql.Tx,
	ctx context.Context,
	id int64,
) error {
	_, err := tx.ExecContext(ctx, srStageDoneSQL, id)
	return err
}

const (
	srStageDoneSQL = `
UPDATE waterway.sounding_results SET staging_done = true
WHERE id = $1`

	insertPointsSQL = `
INSERT INTO waterway.sounding_results (
  bottleneck_id,
  date_info,
  depth_reference,
  point_cloud,
  area
) VALUES (
  (SELECT bottleneck_id from waterway.bottlenecks where objnam = $1),
  $2::date,
  $3,
  ST_Transform(ST_GeomFromWKB($4, $6::integer), 4326)::geography,
  (SELECT CASE WHEN $5::bytea IS NULL THEN
    ST_Transform(
      ST_ConcaveHull(
          ST_Force2D(ST_GeomFromWKB($4, $6::integer)), 0.7), 4326)::geography
  ELSE
    ST_Transform(ST_GeomFromWKB($5, $6::integer), 4326)::geography
  END)
)
RETURNING
  id,
  CASE WHEN ST_Y(ST_Centroid(point_cloud::geometry)) > 0 THEN
    32600
  ELSE
    32700
  END + floor((ST_X(ST_Centroid(point_cloud::geometry))+180)/6)::int + 1`

	insertOctreeSQL = `
UPDATE waterway.sounding_results SET
  octree_checksum = $2, octree_index = $3
WHERE id = $1`

	insertContourSQL = `
INSERT INTO waterway.sounding_results_contour_lines (
  sounding_result_id,
  height,
  lines
)
SELECT
  $1,
  $2,
  ST_Transform(
    ST_Multi(
      ST_CollectionExtract(
        ST_SimplifyPreserveTopology(
          ST_GeomFromWKB($4, $3::integer),
          $5
        ),
        2
      )
    ),
    4326
  )
FROM waterway.sounding_results sr
WHERE id = $1
`
)

func (sr *SoundingResult) FromString(data string) error {
	return json.NewDecoder(strings.NewReader(data)).Decode(sr)
}

func (sr *SoundingResult) ToString() (string, error) {
	var b strings.Builder
	if err := json.NewEncoder(&b).Encode(sr); err != nil {
		return "", err
	}
	return b.String(), nil
}

func (sr *SoundingResult) Do(
	importID int64,
	ctx context.Context,
	conn *sql.Conn,
	feedback Feedback,
) error {

	z, err := zip.OpenReader(filepath.Join(sr.Dir, "sr.zip"))
	if err != nil {
		return err
	}
	defer z.Close()

	feedback.Info("Looking for 'meta.json'")
	mf := common.FindInZIP(z, "meta.json")
	if mf == nil && !sr.completeOverride() {
		return errors.New("Cannot find 'meta.json'")
	}

	m, err := sr.loadMeta(mf)
	if err != nil {
		return err
	}

	if err := m.Validate(conn, ctx); err != nil {
		return common.ToError(err)
	}

	feedback.Info("Looking for '*.xyz'")
	xyzf := common.FindInZIP(z, ".xyz")
	if xyzf == nil {
		return errors.New("Cannot find any *.xyz file")
	}

	xyz, err := loadXYZ(xyzf, feedback)
	if err != nil {
		return err
	}

	if len(xyz) == 0 {
		return errors.New("XYZ does not contain any vertices.")
	}

	// Is there a boundary shapefile in the ZIP archive?
	polygon, err := loadBoundary(z)
	if err != nil {
		return err
	}

	tx, err := conn.BeginTx(ctx, nil)
	if err != nil {
		return err
	}
	defer tx.Rollback()

	var id int64
	var epsg uint32
	start := time.Now()

	err = tx.QueryRow(insertPointsSQL,
		m.Bottleneck,
		m.Date.Time,
		m.DepthReference,
		xyz.AsWKB(),
		polygon.AsWBK(),
		m.EPSG,
	).Scan(&id, &epsg)
	xyz, polygon = nil, nil // not need from now on.
	feedback.Info("storing points took %s", time.Since(start))
	if err != nil {
		return err
	}

	feedback.Info("Best suited UTM EPSG: %d", epsg)

	feedback.Info("Triangulate...")
	start = time.Now()
	tin, err := octree.GenerateTinByID(conn, ctx, id, epsg)
	feedback.Info("triangulation took %s", time.Since(start))
	if err != nil {
		return err
	}

	if tin == nil {
		return errors.New("cannot load TIN from database")
	}

	feedback.Info("Building octree...")
	builder := octree.NewBuilder(tin)
	start = time.Now()
	builder.Build()
	octreeIndex, err := builder.Bytes()
	tin = nil // not needed from now on
	feedback.Info("building octree took %s", time.Since(start))
	if err != nil {
		return err
	}

	feedback.Info("Store octree...")
	start = time.Now()
	h := sha1.New()
	h.Write(octreeIndex)
	checksum := hex.EncodeToString(h.Sum(nil))
	_, err = tx.Exec(insertOctreeSQL, id, checksum, octreeIndex)
	feedback.Info("storing octree index took %s", time.Since(start))
	if err != nil {
		return err
	}

	tree := builder.Tree()
	builder = nil // not needed from now on

	feedback.Info("Generate contour lines...")
	start = time.Now()
	err = generateContours(tree, tx, id)
	feedback.Info("generating and storing contour lines took %s", time.Since(start))
	if err != nil {
		return err
	}

	// Store for potential later removal.
	if err = track(tx, ctx, importID, "waterway.sounding_results", id); err != nil {
		return err
	}

	if err = tx.Commit(); err == nil {
		feedback.Info("Storing sounding result was successful.")
	}
	return err
}

func (sr *SoundingResult) CleanUp() error {
	return os.RemoveAll(sr.Dir)
}

func (sr *SoundingResult) completeOverride() bool {
	return sr.EPSG != nil && sr.Bottleneck != nil && sr.Date != nil && sr.DepthReference != nil
}

func (sr *SoundingResult) loadMeta(f *zip.File) (*models.SoundingResultMeta, error) {
	if f == nil {
		return &models.SoundingResultMeta{
			Date:           *sr.Date,
			Bottleneck:     *sr.Bottleneck,
			EPSG:           *sr.EPSG,
			DepthReference: *sr.DepthReference,
		}, nil
	}
	r, err := f.Open()
	if err != nil {
		return nil, err
	}
	defer r.Close()
	var m models.SoundingResultMeta
	if err := m.Decode(r); err != nil {
		return nil, err
	}

	// Apply overrides
	if sr.Date != nil {
		m.Date = *sr.Date
	}
	if sr.Bottleneck != nil {
		m.Bottleneck = *sr.Bottleneck
	}
	if sr.EPSG != nil {
		m.EPSG = *sr.EPSG
	}
	if sr.DepthReference != nil {
		m.DepthReference = *sr.DepthReference
	}

	return &m, nil
}

func loadXYZReader(r io.Reader, feedback Feedback) (octree.MultiPointZ, error) {
	mpz := make(octree.MultiPointZ, 0, 250000)
	s := bufio.NewScanner(r)

	for line := 1; s.Scan(); line++ {
		text := s.Text()
		var p octree.Vertex
		// fmt.Sscanf(text, "%f,%f,%f") is 4 times slower.
		idx := strings.IndexByte(text, ',')
		if idx == -1 {
			feedback.Warn("format error in line %d", line)
			continue
		}
		var err error
		if p.X, err = strconv.ParseFloat(text[:idx], 64); err != nil {
			feedback.Warn("format error in line %d: %v", line, err)
			continue
		}
		text = text[idx+1:]
		if idx = strings.IndexByte(text, ','); idx == -1 {
			feedback.Warn("format error in line %d", line)
			continue
		}
		if p.Y, err = strconv.ParseFloat(text[:idx], 64); err != nil {
			feedback.Warn("format error in line %d: %v", line, err)
			continue
		}
		text = text[idx+1:]
		if p.Z, err = strconv.ParseFloat(text, 64); err != nil {
			feedback.Warn("format error in line %d: %v", line, err)
			continue
		}
		mpz = append(mpz, p)
	}

	if err := s.Err(); err != nil {
		return nil, err
	}

	return mpz, nil
}

func loadXYZ(f *zip.File, feedback Feedback) (octree.MultiPointZ, error) {
	r, err := f.Open()
	if err != nil {
		return nil, err
	}
	defer r.Close()
	return loadXYZReader(r, feedback)
}

func loadBoundary(z *zip.ReadCloser) (Polygon, error) {
	shpF := common.FindInZIP(z, ".shp")
	if shpF == nil {
		return nil, nil
	}
	prefix := strings.TrimSuffix(shpF.Name, path.Ext(shpF.Name))
	dbfF := common.FindInZIP(z, prefix+".dbf")
	if dbfF == nil {
		return nil, fmt.Errorf("No DBF file found for %s", shpF.Name)
	}

	shpR, err := shpF.Open()
	if err != nil {
		return nil, err
	}

	dbfR, err := dbfF.Open()
	if err != nil {
		shpR.Close()
		return nil, err
	}
	sr := shp.SequentialReaderFromExt(shpR, dbfR)
	defer sr.Close()

	if !sr.Next() {
		return nil, sr.Err()
	}

	_, s := sr.Shape()
	if s == nil {
		return nil, sr.Err()
	}

	return shapeToPolygon(s)
}

func generateContours(tree *octree.Tree, tx *sql.Tx, id int64) error {
	stmt, err := tx.Prepare(insertContourSQL)
	if err != nil {
		return err
	}
	defer stmt.Close()

	// Adjust contour lines heights to multiples of contourStepWidth
	heights := make([]float64, 0)
	h := contourStepWidth * math.Ceil(tree.Min.Z/contourStepWidth)
	for ; h <= tree.Max.Z; h += contourStepWidth {
		heights = append(heights, h)
	}

	octree.DoContours(tree, heights, func(res *octree.ContourResult) {
		if err == nil {
			_, err = stmt.Exec(
				id, res.Height, tree.EPSG,
				res.Lines.AsWKB2D(),
				contourTolerance)
		}
	})

	return err
}