# HG changeset patch
# User Sascha L. Teichmann <sascha.teichmann@intevation.de>
# Date 1551365732 -3600
# Node ID c85b16db8a02caf562e727d11dabc0ee8014935f
# Parent  930ca9c4e2a794de2624d9bb18aa460fcbfe353e
Calculate better triangulation and store it into database.
A bit fake by now. Needs extra:

CREATE TABLE diff_contour_lines (
  height numeric NOT NULL,
  lines geography(MultiLineString,4326) NOT NULL
);

GRANT ALL ON diff_contour_lines TO sophie;

diff -r 930ca9c4e2a7 -r c85b16db8a02 cmd/octreediff/main.go
--- a/cmd/octreediff/main.go	Thu Feb 28 12:34:47 2019 +0100
+++ b/cmd/octreediff/main.go	Thu Feb 28 15:55:32 2019 +0100
@@ -38,6 +38,8 @@
 	second     = flag.String("second", "", "date of the second sounding result")
 )
 
+const contourTolerance = 0.1
+
 const (
 	loadOctreeSQL = `
 SELECT sr.octree_index
@@ -47,8 +49,56 @@
   AND sr.octree_index IS NOT NULL`
 
 	triangulateSQL = `
-SELECT ST_AsBinary(ST_DelaunayTriangles(ST_GeomFromWKB($1, $2::int), 0, 2))
+WITH joined AS (
+  SELECT
+    sr.area    AS area,
+    sr.date_info AS date_info
+  FROM waterway.sounding_results sr JOIN
+     waterway.bottlenecks bn ON sr.bottleneck_id = bn.id
+  WHERE bn.bottleneck_id = $1
+),
+inter AS (
+  SELECT
+  ST_Buffer(
+    ST_intersection(
+      (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $3::date),
+      (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $4::date)
+    ),
+    0.001) AS area
+),
+triangles AS (
+  SELECT t.geom AS geom, ST_MakePolygon(ST_ExteriorRing(t.geom)) AS poly FROM (
+    SELECT (ST_Dump(
+      ST_DelaunayTriangles(ST_GeomFromWKB($5, $2::int), 0, 2))).geom
+    ) t
+)
+SELECT ST_AsBinary(ST_Collect(triangles.geom)) FROM triangles, inter
+WHERE ST_Covers(inter.area, triangles.poly)
 `
+
+	//INSERT INTO redis_diff_countour_lines (
+	//INSERT INTO diff_contour_lines (
+	//INSERT INTO diff_contour_lines_extern (
+	insertContourSQL = `
+INSERT INTO diff_contour_lines (
+  height,
+  lines
+)
+SELECT
+  $1,
+  ST_Transform(
+    ST_Multi(
+      ST_CollectionExtract(
+        ST_SimplifyPreserveTopology(
+          ST_Multi(ST_Collectionextract(
+            ST_MakeValid(ST_GeomFromWKB($2, $3::integer)), 2)),
+          $4
+        ),
+        2
+      )
+    ),
+    4326
+  )`
 )
 
 func check(err error) {
@@ -311,8 +361,11 @@
 		if err := conn.QueryRowContext(
 			ctx,
 			triangulateSQL,
+			bottleneck,
+			first.EPSG,
+			firstDate,
+			secondDate,
 			data,
-			first.EPSG,
 		).Scan(&tin); err != nil {
 			return err
 		}
@@ -361,12 +414,31 @@
 
 		var dataSize int
 
+		stmt, err := conn.PrepareContext(ctx, insertContourSQL)
+		if err != nil {
+			return err
+		}
+		defer stmt.Close()
+
 		octree.DoContours(tree, heights, func(res *octree.ContourResult) {
-			// TODO: Store them.
-			log.Printf("%f: segments: %d\n", res.Height, len(res.Lines))
-			dataSize += len(res.Lines.AsWKB2D())
+			if err == nil && len(res.Lines) > 0 {
+				log.Printf("%f: lines: %d\n", res.Height, len(res.Lines))
+				wkb := res.Lines.AsWKB2D()
+				dataSize += len(wkb)
+				_, err = stmt.ExecContext(
+					ctx,
+					res.Height,
+					wkb,
+					first.EPSG,
+					contourTolerance,
+				)
+			}
 		})
 
+		if err != nil {
+			return err
+		}
+
 		now = time.Now()
 		log.Printf("Number of iso lines: %d\n", len(heights))
 		log.Printf("Total WKB size: %.2fMB\n", float64(dataSize)/(1024*1024))
diff -r 930ca9c4e2a7 -r c85b16db8a02 pkg/octree/contours.go
--- a/pkg/octree/contours.go	Thu Feb 28 12:34:47 2019 +0100
+++ b/pkg/octree/contours.go	Thu Feb 28 15:55:32 2019 +0100
@@ -15,6 +15,7 @@
 package octree
 
 import (
+	"log"
 	"runtime"
 	"sync"
 )
@@ -66,6 +67,12 @@
 // function in order of the original heights values.
 func DoContours(tree *Tree, heights []float64, store func(*ContourResult)) {
 
+	if len(heights) == 0 {
+		return
+	}
+
+	log.Printf("num heights: %d\n", len(heights))
+
 	contours := make([]*ContourResult, len(heights))
 
 	for i, h := range heights {
diff -r 930ca9c4e2a7 -r c85b16db8a02 pkg/octree/tree.go
--- a/pkg/octree/tree.go	Thu Feb 28 12:34:47 2019 +0100
+++ b/pkg/octree/tree.go	Thu Feb 28 15:55:32 2019 +0100
@@ -238,12 +238,14 @@
 			} else {
 				max = mid
 			}
-			if index := ot.index[pos:]; len(index) > 3 {
-				stack = append(stack,
-					frame{index[0], min, max},
-					frame{index[1], min, max},
-					frame{index[2], min, max},
-					frame{index[3], min, max})
+			if pos < int32(len(ot.index)) {
+				if index := ot.index[pos:]; len(index) > 3 {
+					stack = append(stack,
+						frame{index[0], min, max},
+						frame{index[1], min, max},
+						frame{index[2], min, max},
+						frame{index[3], min, max})
+				}
 			}
 		} else { // leaf
 			pos = -pos - 1