Mercurial > gemma
view pkg/octree/areas.go @ 4577:82029885c11b iso-areas
Fixed systematic offset of new iso areas.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Tue, 08 Oct 2019 10:37:02 +0200 |
parents | 6415368c2c60 |
children | fb09a43b062e |
line wrap: on
line source
// 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() quat := 0.25 * cellSize 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 y1 := py - quat y2 := py + quat for j := range row { var n int var sum float64 if v, ok := tree.Value(px-quat, y1); ok { sum = v n = 1 } if v, ok := tree.Value(px-quat, y2); ok { sum += v n++ } if v, ok := tree.Value(px+quat, y1); ok { sum += v n++ } if v, ok := tree.Value(px+quat, y2); ok { sum += v n++ } if n > 0 { row[j] = sum / float64(n) } 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(0.5, min.X, 1.5, min.X+cellSize) reprojY := common.Linear(0.5, min.Y, 1.5, min.Y+cellSize) 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 }