Mercurial > gemma
changeset 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 | f7b57136c800 |
children | a4042ab02e05 |
files | pkg/imports/sr.go |
diffstat | 1 files changed, 223 insertions(+), 12 deletions(-) [+] |
line wrap: on
line diff
--- a/pkg/imports/sr.go Wed Oct 02 16:43:34 2019 +0200 +++ b/pkg/imports/sr.go Wed Oct 02 18:35:09 2019 +0200 @@ -28,15 +28,19 @@ "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 @@ -70,9 +74,17 @@ ) 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" @@ -177,6 +189,31 @@ 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 @@ -619,15 +656,10 @@ return nil, err } feedback.Info("Storing octree index took %s.", time.Since(start)) - feedback.Info("Generate contour lines") - - start = time.Now() - err = generateContours(ctx, tx, feedback, builder.Tree(), id) + err = generateIsos(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 { @@ -827,23 +859,19 @@ 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") @@ -871,6 +899,189 @@ 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(