diff pkg/imports/sr.go @ 2529:45d51a49f191

SR import: Use own triangulation and clipping when importing sounding results.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Wed, 06 Mar 2019 17:51:58 +0100
parents 1ec4c5633eb6
children 3f61b84ae7a6
line wrap: on
line diff
--- a/pkg/imports/sr.go	Wed Mar 06 16:26:45 2019 +0100
+++ b/pkg/imports/sr.go	Wed Mar 06 17:51:58 2019 +0100
@@ -104,36 +104,33 @@
   WHERE import_id = $1 AND
         relation = 'waterway.sounding_results'::regclass)`
 
-	insertPointsSQL = `
+	insertHullSQL = `
 INSERT INTO waterway.sounding_results (
   bottleneck_id,
   date_info,
   depth_reference,
-  point_cloud,
   area
-) VALUES (
+) SELECT
   (SELECT 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)
-)
+  (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,
-  ST_X(ST_Centroid(point_cloud::geometry)),
-  ST_Y(ST_Centroid(point_cloud::geometry)),
-  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`
+  ST_X(ST_Centroid(area::geometry)),
+  ST_Y(ST_Centroid(area::geometry)),
+  best_utm(area::geometry),
+  ST_AsBinary(ST_Transform(area::geometry, best_utm(area::geometry)))
+`
 
+	reprojectPointsSQL = `
+SELECT ST_AsBinary(ST_Transform(ST_GeomFromWKB($1, $2::integer), $3::integer))
+`
 	insertOctreeSQL = `
 UPDATE waterway.sounding_results SET
   octree_checksum = $2, octree_index = $3
@@ -174,6 +171,8 @@
 	feedback Feedback,
 ) (interface{}, error) {
 
+	begin := time.Now()
+
 	z, err := zip.OpenReader(filepath.Join(sr.Dir, "sr.zip"))
 	if err != nil {
 		return nil, err
@@ -234,63 +233,116 @@
 	)
 	start := time.Now()
 
-	err = tx.QueryRow(insertPointsSQL,
+	var hull []byte
+
+	xyzWKB := xyz.AsWKB()
+
+	err = tx.QueryRow(insertHullSQL,
 		m.Bottleneck,
 		m.Date.Time,
 		m.DepthReference,
-		xyz.AsWKB(),
+		xyzWKB,
 		polygon.asWKB(),
 		m.EPSG,
-	).Scan(&id, &lat, &lon, &epsg)
+	).Scan(
+		&id,
+		&lat,
+		&lon,
+		&epsg,
+		&hull,
+	)
 	xyz, polygon = nil, nil // not need from now on.
-	feedback.Info("storing points took %s", time.Since(start))
+	feedback.Info("Calculating hull took %s.", time.Since(start))
+	if err != nil {
+		return nil, err
+	}
+	feedback.Info("Best suited UTM EPSG: %d", epsg)
+
+	start = time.Now()
+
+	var clippingPolygon octree.Polygon
+	if err := clippingPolygon.FromWKB(hull); err != nil {
+		return nil, err
+	}
+	clippingPolygon.Indexify()
+
+	feedback.Info("Building clipping polygon took %v.", time.Since(start))
+
+	start = time.Now()
+
+	var reproj []byte
+
+	err = tx.QueryRow(reprojectPointsSQL,
+		xyzWKB,
+		m.EPSG,
+		epsg,
+	).Scan(&reproj)
+
 	if err != nil {
 		return nil, err
 	}
 
-	feedback.Info("Best suited UTM EPSG: %d", epsg)
+	if err := xyz.FromWKB(reproj); err != nil {
+		return nil, err
+	}
 
-	feedback.Info("Triangulate...")
+	feedback.Info("Reprojecting points took %v.", time.Since(start))
+	feedback.Info("Number of reprojected points: %d", len(xyz))
+
 	start = time.Now()
-	tin, err := octree.GenerateTinByID(ctx, conn, id, epsg)
-	feedback.Info("triangulation took %s", time.Since(start))
+
+	tri, err := octree.Triangulate(xyz)
 	if err != nil {
 		return nil, err
 	}
 
-	if tin == nil {
-		return nil, errors.New("cannot load TIN from database")
-	}
+	feedback.Info("Triangulation took %v.", time.Since(start))
+
+	start = time.Now()
+
+	tin := tri.Tin()
+
+	var str octree.STRTree
+	str.Build(tin)
 
-	feedback.Info("Building octree...")
+	feedback.Info("Building STR tree took %v", time.Since(start))
+
+	start = time.Now()
+
+	removed := str.Clip(&clippingPolygon)
+	feedback.Info("Clipping STR tree took %v.", time.Since(start))
+	feedback.Info("Number of triangles to clip %d.", len(removed))
+
+	start = time.Now()
+
 	builder := octree.NewBuilder(tin)
-	start = time.Now()
-	builder.Build(nil)
+	builder.Build(removed)
 	octreeIndex, err := builder.Bytes()
-	tin = nil // not needed from now on
-	feedback.Info("building octree took %s", time.Since(start))
+	feedback.Info("Building octree took %v.",
+		time.Since(start))
 	if err != nil {
 		return nil, 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))
+	feedback.Info("Storing octree index took %s.",
+		time.Since(start))
 	if err != nil {
 		return nil, err
 	}
 
 	tree := builder.Tree()
+	tree.EPSG = epsg
 	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))
+	feedback.Info("Generating and storing contour lines took %s.",
+		time.Since(start))
 	if err != nil {
 		return nil, err
 	}
@@ -300,10 +352,17 @@
 		return nil, err
 	}
 
-	if err = tx.Commit(); err == nil {
-		feedback.Info("Storing sounding result was successful.")
+	if err = tx.Commit(); err != nil {
+		feedback.Error(
+			"Storing sounding result failed after %v.",
+			time.Since(begin))
+		return nil, err
 	}
 
+	feedback.Info(
+		"Storing sounding result was successful after %v.",
+		time.Since(begin))
+
 	summary := struct {
 		Bottleneck string      `json:"bottleneck"`
 		Date       models.Date `json:"date"`
@@ -473,7 +532,7 @@
 	}
 
 	octree.DoContours(tree, heights, func(res *octree.ContourResult) {
-		if err == nil {
+		if err == nil && len(res.Lines) > 0 {
 			_, err = stmt.Exec(
 				id, res.Height, tree.EPSG,
 				res.Lines.AsWKB2D(),