Mercurial > gemma

changeset 801:b6a1779ffb42

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Removed old cross section controller.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Thu, 27 Sep 2018 13:18:50 +0200
parents 0689f4b7c99f
children 327aa4a18a1c
files pkg/controllers/cross.go pkg/controllers/octreecross.go pkg/controllers/routes.go pkg/models/cross.go
diffstat 4 files changed, 104 insertions(+), 397 deletions(-) [+]
line wrap: on
line diff
--- a/pkg/controllers/cross.go	Thu Sep 27 12:49:45 2018 +0200
+++ b/pkg/controllers/cross.go	Thu Sep 27 13:18:50 2018 +0200
@@ -1,86 +1,141 @@
 package controllers
 
 import (
+	"context"
 	"database/sql"
 	"log"
 	"net/http"
 	"time"
 
 	"gemma.intevation.de/gemma/pkg/models"
+	"gemma.intevation.de/gemma/pkg/octree"
 )
 
-const crossSQL = `
-WITH line AS (
-SELECT ST_LineMerge(ST_Union(ST_3DIntersection(
-  ST_Translate(
-    ST_Extrude(
-      ST_GeomFromWKB($1, 4326),
-    0, 0, 1000),
-   0, 0, -500),
-   geom))) AS geom
-FROM waterway.meshes m JOIN waterway.sounding_results sr ON m.sounding_result_id = sr.id
-WHERE ST_Intersects(geom, ST_GeomFromWKB($1, 4326)) AND
-  sr.bottleneck_id = $2 AND sr.date_info = $3
-)
-SELECT ST_AsBinary((ST_Dump(ST_Intersection(line.geom, sr.area::geometry))).geom)
-  FROM line, waterway.sounding_results sr
-  WHERE sr.bottleneck_id = $2 AND sr.date_info = $3
-`
+const WGS84 = 4326
+
+func reproject(
+	rp *models.Reprojector,
+	src models.GeoJSONLineCoordinates,
+	ctx context.Context,
+) (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(
+			s.Lat, s.Lon, ctx,
+		); err != nil {
+			return nil, err
+		}
+	}
+	return dst, nil
+}
+
+func projectBack(
+	rp *models.Reprojector,
+	lines octree.MultiLineStringZ,
+	ctx context.Context,
+) (models.GeoJSONMultiLineCoordinatesZ, error) {
+
+	out := make(models.GeoJSONMultiLineCoordinatesZ, len(lines))
+
+	for i, segment := range lines {
+
+		coords := make(models.GeoJSONLineCoordinatesZ, len(segment))
+		out[i] = coords
+		for j, v := range segment {
+			lat, lon, err := rp.Reproject(v.X, v.Y, ctx)
+			if err != nil {
+				return nil, err
+			}
+			coords[j] = models.GeoJSONCoordinateZ{Lat: lat, Lon: lon, Z: v.Z}
+		}
+	}
+	return out, nil
+}
 
 func crossSection(
 	input interface{},
 	req *http.Request,
-	db *sql.Conn,
+	conn *sql.Conn,
 ) (jr JSONResult, err error) {
 
 	csi := input.(*models.CrossSectionInput)
 
 	start := time.Now()
 
-	var rows *sql.Rows
-	if rows, err = db.QueryContext(
-		req.Context(),
-		crossSQL,
-		csi.Geometry.Coordinates.AsWKB(),
-		csi.Properties.Bottleneck,
-		csi.Properties.Date.Time,
+	tree, err := octree.Cache.Get(
+		csi.Properties.Bottleneck, csi.Properties.Date.Time,
+		conn, req.Context())
+
+	log.Printf("loading octree took: %s\n", time.Since(start))
+	if err != nil {
+		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(
+		conn, req.Context(),
+		WGS84, tree.EPSG,
 	); err != nil {
 		return
 	}
-	defer rows.Close()
-
-	var segments models.GeoJSONMultiLineCoordinatesZ
+	defer rp.Close()
 
-	for rows.Next() {
-		var segment models.GeoJSONLineCoordinatesZ
-		if err = rows.Scan(&segment); err != nil {
-			return
-		}
-		segments = append(segments, segment)
-	}
+	coords, err := reproject(rp, csi.Geometry.Coordinates, req.Context())
 
-	if err = rows.Err(); err != nil {
+	log.Printf("transforming input coords took: %s\n", time.Since(start))
+	if err != nil {
 		return
 	}
 
-	log.Printf("query took: %v\n", time.Since(start))
-
 	start = time.Now()
 
-	joined := segments.Join()
+	var segments octree.MultiLineStringZ
+
+	for i := 0; i < len(coords)-1; i++ {
+		c1 := &coords[i]
+		c2 := &coords[i+1]
 
-	log.Printf("joining took: %v\n", time.Since(start))
-	log.Printf("segments before/after: %d %d\n", len(segments), len(joined))
+		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 err2 := dumpProfile(
-			csi.Geometry.Coordinates[0],
-			csi.Geometry.Coordinates[len(csi.Geometry.Coordinates)-1],
-			joined,
-		); err2 != nil {
-			log.Printf("error: %v\n", err2)
+		if len(line) > 0 {
+			log.Printf("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("joined length: %d\n", len(joined))
+			segments = append(segments, joined...)
 		}
-	*/
+
+	}
+	log.Printf("octree traversal took: %s\n", time.Since(start))
+
+	// The result have to be WGS84. So project the result back.
+	start = time.Now()
+
+	rp.FromEPSG, rp.ToEPSG = tree.EPSG, WGS84
+
+	var joined models.GeoJSONMultiLineCoordinatesZ
+	joined, err = projectBack(rp, segments, req.Context())
+
+	log.Printf("projecting back took: %s\n", time.Since(start))
+	if err != nil {
+		return
+	}
 
 	jr = JSONResult{
 		Result: &models.CrossSectionOutput{
@@ -95,75 +150,3 @@
 
 	return
 }
-
-/*
-func dumpProfile(
-	start models.GeoJSONCoordinate,
-	stop models.GeoJSONCoordinate,
-	segments models.GeoJSONMultiLineCoordinatesZ,
-) error {
-
-	w, err := os.Create("/tmp/cross.txt")
-	if err != nil {
-		return err
-	}
-	defer w.Close()
-
-	begin := models.GeoJSONCoordinateZ{
-		Lat: start.Lat,
-		Lon: start.Lon,
-	}
-	end := models.GeoJSONCoordinateZ{
-		Lat: stop.Lat,
-		Lon: stop.Lon,
-	}
-
-	last := begin
-
-	var total float64
-
-	minz := 10000.0
-
-	for _, line := range segments {
-		for _, coord := range line {
-			if coord.Z < minz {
-				minz = coord.Z
-			}
-			total += last.Distance(coord)
-			last = coord
-		}
-	}
-
-	log.Printf("length: %.3f minz:  %f\n", total, minz)
-
-	var pos float64
-
-	for i, line := range segments {
-		for j, coord := range line {
-			if i == 0 && j == 0 {
-				if !begin.Equals(coord) {
-					pos = begin.Distance(coord)
-					fmt.Fprintf(w, "%.3f %f\n", 0.0, 200.0)
-					fmt.Fprintf(w, "%.3f %f\n", pos-0.0001, 200.0)
-					fmt.Fprintf(w, "%.3f %f\n", pos, coord.Z-minz)
-					continue
-				}
-			} else if j == 0 {
-				fmt.Fprintf(w, "%.3f %f\n", pos+0.0001, 200.0)
-				fmt.Fprintf(w, "%.3f %f\n", pos+last.Distance(coord)-0.0001, 200.0)
-				continue
-			}
-			pos += last.Distance(coord)
-			fmt.Fprintf(w, "%.3f %f\n", pos, coord.Z-minz)
-			last = coord
-		}
-	}
-
-	if !last.Equals(end) {
-		fmt.Fprintf(w, "%.3f %f\n", pos+0.0001, 200.0)
-		fmt.Fprintf(w, "%.3f %f\n", pos+last.Distance(end), 200.0)
-	}
-
-	return nil
-}
-*/
--- a/pkg/controllers/octreecross.go	Thu Sep 27 12:49:45 2018 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,152 +0,0 @@
-package controllers
-
-import (
-	"context"
-	"database/sql"
-	"log"
-	"net/http"
-	"time"
-
-	"gemma.intevation.de/gemma/pkg/models"
-	"gemma.intevation.de/gemma/pkg/octree"
-)
-
-const WGS84 = 4326
-
-func reproject(
-	rp *models.Reprojector,
-	src models.GeoJSONLineCoordinates,
-	ctx context.Context,
-) (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(
-			s.Lat, s.Lon, ctx,
-		); err != nil {
-			return nil, err
-		}
-	}
-	return dst, nil
-}
-
-func projectBack(
-	rp *models.Reprojector,
-	lines octree.MultiLineStringZ,
-	ctx context.Context,
-) (models.GeoJSONMultiLineCoordinatesZ, error) {
-
-	out := make(models.GeoJSONMultiLineCoordinatesZ, len(lines))
-
-	for i, segment := range lines {
-
-		coords := make(models.GeoJSONLineCoordinatesZ, len(segment))
-		out[i] = coords
-		for j, v := range segment {
-			lat, lon, err := rp.Reproject(v.X, v.Y, ctx)
-			if err != nil {
-				return nil, err
-			}
-			coords[j] = models.GeoJSONCoordinateZ{Lat: lat, Lon: lon, Z: v.Z}
-		}
-	}
-	return out, nil
-}
-
-func octreeCrossSection(
-	input interface{},
-	req *http.Request,
-	conn *sql.Conn,
-) (jr JSONResult, err error) {
-
-	csi := input.(*models.CrossSectionInput)
-
-	start := time.Now()
-
-	tree, err := octree.Cache.Get(
-		csi.Properties.Bottleneck, csi.Properties.Date.Time,
-		conn, req.Context())
-
-	log.Printf("loading octree took: %s\n", time.Since(start))
-	if err != nil {
-		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(
-		conn, req.Context(),
-		WGS84, tree.EPSG,
-	); err != nil {
-		return
-	}
-	defer rp.Close()
-
-	coords, err := reproject(rp, csi.Geometry.Coordinates, req.Context())
-
-	log.Printf("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("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("joined length: %d\n", len(joined))
-			segments = append(segments, joined...)
-		}
-
-	}
-	log.Printf("octree traversal took: %s\n", time.Since(start))
-
-	// The result have to be WGS84. So project the result back.
-	start = time.Now()
-
-	rp.FromEPSG, rp.ToEPSG = tree.EPSG, WGS84
-
-	var joined models.GeoJSONMultiLineCoordinatesZ
-	joined, err = projectBack(rp, segments, req.Context())
-
-	log.Printf("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{}{},
-		},
-	}
-
-	return
-}
--- a/pkg/controllers/routes.go	Thu Sep 27 12:49:45 2018 +0200
+++ b/pkg/controllers/routes.go	Thu Sep 27 13:18:50 2018 +0200
@@ -105,11 +105,6 @@
 
 	// Cross sections
 
-	api.Handle("/octcross", any(&JSONHandler{
-		Input:  func() interface{} { return new(models.CrossSectionInput) },
-		Handle: octreeCrossSection,
-	})).Methods(http.MethodPost)
-
 	api.Handle("/cross", any(&JSONHandler{
 		Input:  func() interface{} { return new(models.CrossSectionInput) },
 		Handle: crossSection,
--- a/pkg/models/cross.go	Thu Sep 27 12:49:45 2018 +0200
+++ b/pkg/models/cross.go	Thu Sep 27 13:18:50 2018 +0200
@@ -7,7 +7,6 @@
 	"errors"
 	"fmt"
 	"math"
-	"sort"
 	"time"
 )
 
@@ -257,121 +256,3 @@
 	x := dLng * math.Cos(deg2rad((cz.Lat+other.Lat)/2.0))
 	return math.Sqrt(dLat*dLat+x*x) * EarthRadius
 }
-
-type point struct {
-	x float64
-	y float64
-}
-
-func (cz GeoJSONCoordinateZ) quant() point {
-	const xyScale = 1e6
-	return point{
-		math.Round(cz.Lon*xyScale) / xyScale,
-		math.Round(cz.Lat*xyScale) / xyScale,
-	}
-}
-
-func (lcz GeoJSONLineCoordinatesZ) join(other GeoJSONLineCoordinatesZ) GeoJSONLineCoordinatesZ {
-	l1, l2 := len(lcz), len(other)
-	n := make(GeoJSONLineCoordinatesZ, l1+l2-1)
-	copy(n, lcz[:l1-1])
-	n[l1-1] = lcz[l1-1].mid(other[0])
-	copy(n[l1:], other[1:])
-	return n
-}
-
-func (ml GeoJSONMultiLineCoordinatesZ) Join() GeoJSONMultiLineCoordinatesZ {
-
-	if len(ml) < 2 {
-		return ml
-	}
-
-	type value struct {
-		coords GeoJSONLineCoordinatesZ
-		order  int
-	}
-
-	heads := make(map[point]*value)
-
-	for i, coords := range ml {
-		key := coords[len(coords)-1].quant()
-		if v := heads[key]; v != nil {
-			delete(heads, key)
-			v.coords = coords.join(v.coords)
-			heads[v.coords[0].quant()] = v
-		} else {
-			heads[coords[0].quant()] = &value{coords, i}
-		}
-	}
-
-again:
-	for k, v1 := range heads {
-		key := v1.coords[len(v1.coords)-1].quant()
-		if k == key { // cycle detected
-			continue
-		}
-		if v2 := heads[key]; v2 != nil {
-			delete(heads, key)
-			v1.coords = v1.coords.join(v2.coords)
-			goto again
-		}
-	}
-
-	// To make it deterministic sort in input order.
-	tmp := make([]*value, len(heads))
-	var i int
-	for _, v := range heads {
-		tmp[i] = v
-		i++
-	}
-	sort.Slice(tmp, func(i, j int) bool { return tmp[i].order < tmp[j].order })
-
-	out := make(GeoJSONMultiLineCoordinatesZ, len(tmp))
-	for i, v := range tmp {
-		out[i] = v.coords
-	}
-
-	return out
-}
-
-/*
-func (ml GeoJSONMultiLineCoordinatesZ) Join() GeoJSONMultiLineCoordinatesZ {
-
-	if len(ml) < 2 {
-		return ml
-	}
-
-	lists := make(GeoJSONMultiLineCoordinatesZ, len(ml))
-	copy(lists, ml)
-
-next:
-	for j := 0; j < len(lists); {
-		front := lists[j]
-		tail := front[len(front)-1]
-
-		for i := 0; i < len(lists); i++ {
-			if i == j {
-				continue
-			}
-			curr := lists[i]
-			head := curr[0]
-			if tail.Equals(head) {
-				n := make(GeoJSONLineCoordinatesZ, len(front)+len(curr)-1)
-				copy(n, front[:len(front)-1])
-				n[len(front)-1] = head.mid(tail)
-				copy(n[len(front):], curr[1:])
-				lists[j] = n
-				if i < len(lists)-1 {
-					copy(lists[i:], lists[i+1:])
-				}
-				lists[len(lists)-1] = nil
-				lists = lists[:len(lists)-1]
-				continue next
-			}
-		}
-		j++
-	}
-
-	return lists
-}
-*/