Mercurial > gemma
diff pkg/octree/areas.go @ 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 | |
children | 6415368c2c60 |
line wrap: on
line diff
--- /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 +}