# HG changeset patch
# User Sascha L. Teichmann <sascha.teichmann@intevation.de>
# Date 1570125113 -7200
# Node ID 6b107d4e6810fb0ef12110f5811d1e2b7b59b18e
# Parent  a4042ab02e05901e0c342a82584a7b4637a442a8
Add more debug output to trace problem why reprojected iso areas are not in expected bounds.

diff -r a4042ab02e05 -r 6b107d4e6810 pkg/imports/sr.go
--- a/pkg/imports/sr.go	Thu Oct 03 10:05:01 2019 +0200
+++ b/pkg/imports/sr.go	Thu Oct 03 19:51:53 2019 +0200
@@ -24,6 +24,7 @@
 	"errors"
 	"fmt"
 	"io"
+	"log"
 	"math"
 	"os"
 	"path"
@@ -206,7 +207,7 @@
       ST_CollectionExtract(
         ST_SimplifyPreserveTopology(
           ST_Multi(ST_Collectionextract(
-            ST_MakeValid(ST_GeomFromWKB($4, $3::integer)), 2)),
+            ST_MakeValid(ST_GeomFromWKB($4, $3::integer)), 3)),
           $5
         ),
         3
@@ -929,12 +930,12 @@
 	width := max.X - min.X
 	height := max.Y - min.Y
 
-	feedback.Info("width/height: %.2f / %.2f", width, height)
+	feedback.Info("Width/Height: %.2f / %.2f", width, height)
 
 	xcells := int(math.Ceil(width / isoCellSize))
 	ycells := int(math.Ceil(height / isoCellSize))
 
-	feedback.Info("raster size: (%d, %d)", xcells, ycells)
+	feedback.Info("Raster size: (%d, %d)", xcells, ycells)
 
 	// Add border for closing
 	raster := make([]float64, (xcells+2)*(ycells+2))
@@ -982,10 +983,9 @@
 
 	wg.Wait()
 
-	feedback.Info("Rasterizing took %v", time.Since(start))
+	feedback.Info("Rasterizing took %v.", time.Since(start))
 
-	// TODO: Implement me!
-	feedback.Info("Trace iso areas")
+	feedback.Info("Trace iso areas.")
 
 	start = time.Now()
 
@@ -1004,23 +1004,54 @@
 		for hIdx := range cnts {
 			c := tracer.Contours(heights[hIdx])
 
+			if len(c) == 0 {
+				continue
+			}
+
 			// We need to bring it back to the
 			// none raster coordinate system.
 			a := make(wkb.MultiPolygonGeom, len(c))
+
+			rXMin, rXMax := math.MaxFloat64, -math.MaxFloat64
+			rYMin, rYMax := math.MaxFloat64, -math.MaxFloat64
+
 			for i, pl := range c {
 				shell := make(wkb.LinearRingGeom, len(pl))
 				for j, pt := range pl {
+					x := reprojX(pt.X)
+					y := reprojY(pt.Y)
+
+					if x < rXMin {
+						rXMin = x
+					}
+					if y < rYMin {
+						rYMin = y
+					}
+					if x > rXMax {
+						rXMax = x
+					}
+					if y > rYMax {
+						rYMax = y
+					}
+
 					shell[j] = wkb.PointGeom{
-						X: reprojX(pt.X),
-						Y: reprojY(pt.Y),
+						X: x,
+						Y: y,
 					}
 				}
 				a[i] = wkb.PolygonGeom{shell}
 			}
+			log.Printf("(%.2f, %.2f) - (%.2f, %.2f)\n",
+				rXMin, rYMin, rXMax, rYMax)
+			log.Printf("%.2f / %.2f\n", rXMax-rXMin, rYMax-rYMin)
+
 			areas[hIdx] = a
 		}
 	}
 
+	log.Printf("[%.2f, %.2f] - [%.2f, %.2f]\n",
+		min.X, min.Y, max.X, max.Y)
+
 	for n := runtime.NumCPU(); n >= 1; n-- {
 		wg.Add(1)
 		go doContours()
@@ -1033,7 +1064,7 @@
 
 	wg.Wait()
 
-	feedback.Info("Tracing areas took %v", time.Since(start))
+	feedback.Info("Tracing areas took %v.", time.Since(start))
 
 	// Raster and tracer are not needed any more.
 	raster, tracer = nil, nil
@@ -1052,13 +1083,15 @@
 	heights []float64,
 	id int64,
 ) error {
-	feedback.Info("Storing iso areas")
+	feedback.Info("Store iso areas.")
 	total := time.Now()
 	defer func() {
-		feedback.Info("Storing iso areas took %v",
+		feedback.Info("Storing iso areas took %v.",
 			time.Since(total))
 	}()
 
+	log.Printf("EPSG: %d\n", epsg)
+
 	stmt, err := tx.PrepareContext(ctx, insertIsoAreasSQL)
 	if err != nil {
 		return err