view pkg/imports/sr.go @ 971:f9fb6c399f3f

Moved generating of tins to octree package.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Thu, 18 Oct 2018 11:52:13 +0200
parents 1e2dce348cfb
children 17a03a84b0e8
line wrap: on
line source

package imports

import (
	"archive/zip"
	"bufio"
	"bytes"
	"context"
	"database/sql"
	"encoding/binary"
	"encoding/json"
	"errors"
	"fmt"
	"io"
	"log"
	"math"
	"os"
	"path"
	"strconv"
	"strings"
	"time"

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

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

type SoundingResult struct {
	who string
	zip string
}

const SoundingResultDateFormat = "2006-01-02"

type (
	SoundingResultDate struct{ time.Time }

	SoundingResultMeta struct {
		Date           SoundingResultDate `json:"date"`
		Bottleneck     string             `json:"bottleneck"`
		EPSG           uint               `json:"epsg"`
		DepthReference string             `json:"depth-reference"`
	}

	Point struct {
		X float64
		Y float64
	}
	LineString []Point
	Polygon    []LineString
)

const (
	wkbNDR     byte   = 1
	wkbPolygon uint32 = 3
)

const (
	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 $5 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`
)

func (srd *SoundingResultDate) UnmarshalJSON(data []byte) error {
	var s string
	if err := json.Unmarshal(data, &s); err != nil {
		return err
	}
	d, err := time.Parse(SoundingResultDateFormat, s)
	if err == nil {
		*srd = SoundingResultDate{d}
	}
	return err
}

func (sr *SoundingResult) Who() string {
	return sr.who
}

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

func find(needle string, haystack []*zip.File) *zip.File {
	needle = strings.ToLower(needle)
	for _, straw := range haystack {
		if strings.HasSuffix(strings.ToLower(straw.Name), needle) {
			return straw
		}
	}
	return nil
}

func loadMeta(f *zip.File) (*SoundingResultMeta, error) {
	r, err := f.Open()
	if err != nil {
		return nil, err
	}
	defer r.Close()
	var m SoundingResultMeta
	err = json.NewDecoder(r).Decode(&m)
	return &m, err
}

func (m *SoundingResultMeta) validate(conn *sql.Conn) error {

	var b bool
	err := conn.QueryRowContext(context.Background(),
		`SELECT true FROM internal.depth_references WHERE depth_reference = $1`,
		m.DepthReference).Scan(&b)
	switch {
	case err == sql.ErrNoRows:
		return fmt.Errorf("Unknown depth reference '%s'\n", m.DepthReference)
	case err != nil:
		return err
	case !b:
		return errors.New("Unexpected depth reference")
	}

	err = conn.QueryRowContext(context.Background(),
		`SELECT true FROM waterway.bottlenecks WHERE objnam = $1`,
		m.Bottleneck).Scan(&b)
	switch {
	case err == sql.ErrNoRows:
		return fmt.Errorf("Unknown bottleneck '%s'\n", m.Bottleneck)
	case err != nil:
		return err
	case !b:
		return errors.New("Unexpected bottleneck")
	}

	return nil
}

func loadXYZReader(r io.Reader) (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 {
			log.Printf("format error in line %d\n", line)
			continue
		}
		var err error
		if p.X, err = strconv.ParseFloat(text[:idx], 64); err != nil {
			log.Printf("format error in line %d: %v\n", line, err)
			continue
		}
		text = text[idx+1:]
		if idx = strings.IndexByte(text, ','); idx == -1 {
			log.Printf("format error in line %d\n", line)
			continue
		}
		if p.Y, err = strconv.ParseFloat(text[:idx], 64); err != nil {
			log.Printf("format error in line %d: %v\n", line, err)
			continue
		}
		text = text[idx+1:]
		if p.Z, err = strconv.ParseFloat(text, 64); err != nil {
			log.Printf("format error in line %d: %v\n", line, err)
			continue
		}
		mpz = append(mpz, p)
	}

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

	return mpz, nil
}

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

func toPolygon(numParts int32, parts []int32, points []shp.Point) Polygon {
	out := make(Polygon, numParts)
	pos := 0
	for i := range out {
		ps := parts[i]
		line := make(LineString, ps)
		for j := int32(0); j < ps; j, pos = j+1, pos+1 {
			p := &points[pos]
			line[j] = Point{p.X, p.Y}
		}
		out[i] = line
	}
	return out
}

func loadBoundary(z *zip.ReadCloser) (Polygon, error) {
	shpF := find(".shp", z.File)
	if shpF == nil {
		return nil, nil
	}
	prefix := strings.TrimSuffix(shpF.Name, path.Ext(shpF.Name))
	dbfF := find(prefix+".dbf", z.File)
	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()
	}

	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 (sr *SoundingResult) Do(conn *sql.Conn) error {

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

	mf := find("meta.json", z.File)
	if mf == nil {
		return errors.New("Cannot find 'meta.json'")
	}

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

	if err := m.validate(conn); err != nil {
		return err
	}

	xyzf := find(".xyz", z.File)
	if xyzf == nil {
		return errors.New("Cannot find any *.xyz file")
	}

	xyz, err := loadXYZ(xyzf)
	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(context.Background(), nil)
	if err != nil {
		return err
	}
	defer tx.Rollback()

	var id int64
	var epsg uint32

	if err := tx.QueryRow(insertPointsSQL,
		m.Bottleneck,
		m.Date.Time,
		m.DepthReference,
		xyz.AsWKB(),
		polygon.AsWBK(),
		m.EPSG,
	).Scan(&id, &epsg); err != nil {
		return err
	}

	log.Printf("EPSG: %d\n", epsg)

	// TODO: Build octree
	// TODO: Generate iso-lines

	return tx.Commit()
}

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

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

	binary.Write(buf, binary.LittleEndian, wkbNDR)
	binary.Write(buf, binary.LittleEndian, wkbPolygon)
	binary.Write(buf, binary.LittleEndian, uint32(len(p)))

	for _, ring := range p {
		binary.Write(buf, binary.LittleEndian, uint32(len(ring)))
		for _, v := range ring {
			binary.Write(buf, binary.LittleEndian, math.Float64bits(v.X))
			binary.Write(buf, binary.LittleEndian, math.Float64bits(v.Y))
		}
	}

	return buf.Bytes()
}