Mercurial > gemma
changeset 4570:4b3a298b94f8 iso-areas
Moved area tracing to octree package.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Sun, 06 Oct 2019 11:12:32 +0200 |
parents | efb7851f54dd |
children | a413b4a89bcc |
files | pkg/imports/sr.go pkg/octree/areas.go |
diffstat | 2 files changed, 147 insertions(+), 121 deletions(-) [+] |
line wrap: on
line diff
--- a/pkg/imports/sr.go Sun Oct 06 10:29:43 2019 +0200 +++ b/pkg/imports/sr.go Sun Oct 06 11:12:32 2019 +0200 @@ -28,13 +28,10 @@ "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" @@ -925,124 +922,7 @@ 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)) - - feedback.Info("Trace iso areas.") - - start = time.Now() - - tracer := contourmap.FromFloat64s(xcells+2, ycells+2, raster) - - areas := make([]wkb.MultiPolygonGeom, len(heights)) - - // TODO: Check if this correct! - 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]) - - if len(c) == 0 { - continue - } - - // We need to bring it back to the - // none raster coordinate system. - a := make(wkb.MultiPolygonGeom, len(c)) - - for i, pl := range c { - shell := make(wkb.LinearRingGeom, len(pl)) - for j, pt := range pl { - shell[j] = wkb.PointGeom{ - X: reprojX(pt.X), - Y: reprojY(pt.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 + areas := tree.TraceAreas(heights, isoCellSize) return storeAreas( ctx, tx, feedback,
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/octree/areas.go Sun Oct 06 11:12:32 2019 +0200 @@ -0,0 +1,146 @@ +// This is Free Software under GNU Affero General Public License v >= 3.0 +// without warranty, see README.md and license for details. +// +// SPDX-License-Identifier: AGPL-3.0-or-later +// License-Filename: LICENSES/AGPL-3.0.txt +// +// Copyright (C) 2018 by via donau +// – Österreichische Wasserstraßen-Gesellschaft mbH +// Software engineering by Intevation GmbH +// +// Author(s): +// * Sascha L. Teichmann <sascha.teichmann@intevation.de> + +package octree + +import ( + "log" + "math" + "runtime" + "sync" + "time" + + "github.com/fogleman/contourmap" + + "gemma.intevation.de/gemma/pkg/common" + "gemma.intevation.de/gemma/pkg/wkb" +) + +func (tree *Tree) TraceAreas( + heights []float64, + cellSize float64, +) []wkb.MultiPolygonGeom { + min, max := tree.Min, tree.Max + + width := max.X - min.X + height := max.Y - min.Y + + log.Printf("info: Width/Height: %.2f / %.2f\n", width, height) + + xcells := int(math.Ceil(width / cellSize)) + ycells := int(math.Ceil(height / cellSize)) + + log.Printf("info: Raster size: (%d, %d)\n", xcells, ycells) + + start := time.Now() + + // 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)*cellSize + cellSize/2 + px := min.X + cellSize/2 + for j := range row { + v, ok := tree.Value(px, py) + if ok { + row[j] = v + } + px += cellSize + } + } + } + + for n := runtime.NumCPU(); n >= 1; n-- { + wg.Add(1) + go rasterRow() + } + + for i := 0; i < ycells; i++ { + rows <- i + } + close(rows) + + wg.Wait() + log.Printf("info: Rastering took %v\n", time.Since(start)) + + start = time.Now() + + tracer := contourmap.FromFloat64s(xcells+2, ycells+2, raster) + + areas := make([]wkb.MultiPolygonGeom, len(heights)) + + // TODO: Check if this correct! + 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]) + + if len(c) == 0 { + continue + } + + // We need to bring it back to the + // none raster coordinate system. + a := make(wkb.MultiPolygonGeom, len(c)) + + for i, pl := range c { + shell := make(wkb.LinearRingGeom, len(pl)) + for j, pt := range pl { + shell[j] = wkb.PointGeom{ + X: reprojX(pt.X), + Y: reprojY(pt.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() + log.Printf("info: Tracing areas took %v\n", time.Since(start)) + + return areas +}