gemma: pkg/imports/sr.go comparison

comparison pkg/imports/sr.go @ 4562:5cc4042cf07c iso-areas

Started with integrating iso area generation into gemma server. WIP

author	Sascha L. Teichmann <sascha.teichmann@intevation.de>
date	Wed, 02 Oct 2019 18:35:09 +0200
parents	4ca884dfc470
children	a4042ab02e05

comparison

equal deleted inserted replaced

-:f7b57136c800
+:5cc4042cf07c
 	"io"
 	"math"
 	"os"
 	"path"
 	"path/filepath"
+	"runtime"
 	"strconv"
 	"strings"
+	"sync"
 	"time"
+	"github.com/fogleman/contourmap"
 	shp "github.com/jonas-p/go-shp"
 	"gemma.intevation.de/gemma/pkg/common"
 	"gemma.intevation.de/gemma/pkg/models"
 	"gemma.intevation.de/gemma/pkg/octree"
+	"gemma.intevation.de/gemma/pkg/wkb"
 )
 // SoundingResult is a Job to import sounding reults
 // from a ZIP file into the database.
 type SoundingResult struct {
 	contourStepWidth = 0.1
 	contourTolerance = 0.1
 )
 const (
+	// pointsPerSquareMeter is the average number of points
+	// when generating a artifical height model for single beam scans.
 	pointsPerSquareMeter = 2
+)
+const (
+	// isoCellSize is the side length of a raster cell when tracing
+	// iso areas.
+	isoCellSize = 0.5
 )
 // SRJobKind is the unique name of this import job type.
 const SRJobKind JobKind = "sr"
 4326
 )
 FROM waterway.sounding_results sr
 WHERE id = $1
 `
+	insertIsoAreasSQL = `
+INSERT INTO waterway.sounding_results_iso_areas (
+sounding_result_id,
+height,
+areas
+)
+SELECT
+$1,
+$2,
+ST_Transform(
+ST_Multi(
+ST_CollectionExtract(
+ST_SimplifyPreserveTopology(
+ST_Multi(ST_Collectionextract(
+ST_MakeValid(ST_GeomFromWKB($4, $3::integer)), 2)),
+$5
+),
+2
+)
+),
+4326
+)
+FROM waterway.sounding_results sr
+WHERE id = $1
+`
 	selectGaugeLDCSQL = `
 SELECT
 grwl.value,
 grwl.depth_reference
 	_, err = tx.ExecContext(ctx, insertOctreeSQL, id, checksum, octreeIndex)
 	if err != nil {
 		return nil, err
 	}
 	feedback.Info("Storing octree index took %s.", time.Since(start))
-	feedback.Info("Generate contour lines")
+	err = generateIsos(ctx, tx, feedback, builder.Tree(), id)
+	if err != nil {
-	start = time.Now()
+		return nil, err
-	err = generateContours(ctx, tx, feedback, builder.Tree(), id)
+	}
-	if err != nil {
-		return nil, err
-	}
-	feedback.Info("Generating contour lines took %s.",
-		time.Since(start))
 	// Store for potential later removal.
 	if err = track(ctx, tx, importID, "waterway.sounding_results", id); err != nil {
 		return nil, err
 	}
 	}
 	return shapeToPolygon(s)
 }
-func generateContours(
+func generateIsos(
 	ctx context.Context,
 	tx *sql.Tx,
 	feedback Feedback,
 	tree *octree.Tree,
 	id int64,
 ) error {
-	stmt, err := tx.PrepareContext(ctx, insertContourSQL)
-	if err != nil {
-		return err
-	}
-	defer stmt.Close()
 	heights, err := octree.LoadClassBreaks(
 		ctx, tx,
 		"morphology_classbreaks",
 	)
 	if err != nil {
 		feedback.Warn("Loading class breaks failed: %v", err)
 		feedback.Info("Using default class breaks")
 		heights = nil
 		h := contourStepWidth * math.Ceil(tree.Min.Z/contourStepWidth)
 		log.Printf("%.2f - %.2f\n", tree.Min.Z, tree.Max.Z)
 	*/
 	heights = common.DedupFloat64s(heights)
+	if err := generateContours(ctx, tx, feedback, tree, heights, id); err != nil {
+		return err
+	}
+	return generateIsoAreas(ctx, tx, feedback, tree, heights, id)
+}
+func generateIsoAreas(
+	ctx context.Context,
+	tx *sql.Tx,
+	feedback Feedback,
+	tree *octree.Tree,
+	heights []float64,
+	id int64,
+) error {
+	feedback.Info("Generate iso areas")
+	total := time.Now()
+	defer func() {
+		feedback.Info("Generating iso areas took %s.",
+			time.Since(total))
+	}()
+	min, max := tree.Min, tree.Max
+	width := max.X - min.X
+	height := max.Y - min.Y
+	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)
+	// Add border for closing
+	raster := make([]float64, (xcells+2)*(ycells+2))
+	// prefill for no data
+	const nodata = -math.MaxFloat64
+	for i := range raster {
+		raster[i] = nodata
+	}
+	// rasterize the height model
+	var wg sync.WaitGroup
+	rows := make(chan int)
+	rasterRow := func() {
+		defer wg.Done()
+		for i := range rows {
+			pos := (i+1)*(xcells+2) + 1
+			row := raster[pos : pos+xcells]
+			py := min.Y + float64(i)*isoCellSize + isoCellSize/2
+			px := min.X + isoCellSize/2
+			for j := range row {
+				v, ok := tree.Value(px, py)
+				if ok {
+					row[j] = v
+				}
+				px += isoCellSize
+			}
+		}
+	}
+	start := time.Now()
+	for n := runtime.NumCPU(); n >= 1; n-- {
+		wg.Add(1)
+		go rasterRow()
+	}
+	for i := 0; i < ycells; i++ {
+		rows <- i
+	}
+	close(rows)
+	wg.Wait()
+	feedback.Info("Rasterizing took %v", time.Since(start))
+	// TODO: Implement me!
+	feedback.Info("Trace iso areas")
+	start = time.Now()
+	tracer := contourmap.FromFloat64s(xcells+2, ycells+2, raster)
+	areas := make([]wkb.MultiPolygonGeom, len(heights))
+	reprojX := common.Linear(1, min.X, float64(xcells+1), max.X)
+	reprojY := common.Linear(1, min.Y, float64(ycells+1), max.Y)
+	cnts := make(chan int)
+	doContours := func() {
+		defer wg.Done()
+		for hIdx := range cnts {
+			c := tracer.Contours(heights[hIdx])
+			a := make(wkb.MultiPolygonGeom, len(c))
+			for i, pl := range c {
+				shell := make(wkb.LinearRingGeom, len(pl))
+				for j, pt := range pl {
+					x := reprojX(pt.X)
+					y := reprojY(pt.Y)
+					shell[j] = wkb.PointGeom{X: x, Y: y}
+				}
+				a[i] = wkb.PolygonGeom{shell}
+			}
+			areas[hIdx] = a
+		}
+	}
+	for n := runtime.NumCPU(); n >= 1; n-- {
+		wg.Add(1)
+		go doContours()
+	}
+	for i := range heights {
+		cnts <- i
+	}
+	close(cnts)
+	wg.Wait()
+	feedback.Info("Tracing areas took %v", time.Since(start))
+	// Raster and tracer are not needed any more.
+	raster, tracer = nil, nil
+	return storeAreas(ctx, tx, feedback, areas, heights)
+}
+func storeAreas(
+	ctx context.Context,
+	tx *sql.Tx,
+	feedback Feedback,
+	areas []wkb.MultiPolygonGeom,
+	heights []float64,
+) error {
+	feedback.Info("Storing iso areas")
+	total := time.Now()
+	defer func() {
+		feedback.Info("Storing iso areas took %v",
+			time.Since(total))
+	}()
+	stmt, err := tx.PrepareContext(ctx, insertIsoAreasSQL)
+	if err != nil {
+		return err
+	}
+	defer stmt.Close()
+	// TODO: Implement me!
+	return nil
+}
+func generateContours(
+	ctx context.Context,
+	tx *sql.Tx,
+	feedback Feedback,
+	tree *octree.Tree,
+	heights []float64,
+	id int64,
+) error {
+	feedback.Info("Generate contour lines")
+	start := time.Now()
+	defer func() {
+		feedback.Info("Generating contour lines took %v",
+			time.Since(start))
+	}()
+	stmt, err := tx.PrepareContext(ctx, insertContourSQL)
+	if err != nil {
+		return err
+	}
+	defer stmt.Close()
 	octree.DoContours(tree, heights, func(res *octree.ContourResult) {
 		if err == nil && len(res.Lines) > 0 {
 			_, err = stmt.ExecContext(
 				ctx,
 				id, res.Height, tree.EPSG,

Mercurial > gemma

comparison pkg/imports/sr.go @ 4562:5cc4042cf07c iso-areas