Mercurial > gemma

comparison pkg/octree/areas.go @ 4570:4b3a298b94f8 iso-areas

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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
comparison
equal deleted inserted replaced
4569:efb7851f54dd 4570:4b3a298b94f8
1 // This is Free Software under GNU Affero General Public License v >= 3.0
2 // without warranty, see README.md and license for details.
3 //
4 // SPDX-License-Identifier: AGPL-3.0-or-later
5 // License-Filename: LICENSES/AGPL-3.0.txt
6 //
7 // Copyright (C) 2018 by via donau
8 // – Österreichische Wasserstraßen-Gesellschaft mbH
9 // Software engineering by Intevation GmbH
10 //
11 // Author(s):
12 // * Sascha L. Teichmann <sascha.teichmann@intevation.de>
13
14 package octree
15
16 import (
17 "log"
18 "math"
19 "runtime"
20 "sync"
21 "time"
22
23 "github.com/fogleman/contourmap"
24
25 "gemma.intevation.de/gemma/pkg/common"
26 "gemma.intevation.de/gemma/pkg/wkb"
27 )
28
29 func (tree *Tree) TraceAreas(
30 heights []float64,
31 cellSize float64,
32 ) []wkb.MultiPolygonGeom {
33 min, max := tree.Min, tree.Max
34
35 width := max.X - min.X
36 height := max.Y - min.Y
37
38 log.Printf("info: Width/Height: %.2f / %.2f\n", width, height)
39
40 xcells := int(math.Ceil(width / cellSize))
41 ycells := int(math.Ceil(height / cellSize))
42
43 log.Printf("info: Raster size: (%d, %d)\n", xcells, ycells)
44
45 start := time.Now()
46
47 // Add border for closing
48 raster := make([]float64, (xcells+2)*(ycells+2))
49
50 // prefill for no data
51 const nodata = -math.MaxFloat64
52 for i := range raster {
53 raster[i] = nodata
54 }
55
56 // rasterize the height model
57
58 var wg sync.WaitGroup
59
60 rows := make(chan int)
61
62 rasterRow := func() {
63 defer wg.Done()
64 for i := range rows {
65 pos := (i+1)*(xcells+2) + 1
66 row := raster[pos : pos+xcells]
67 py := min.Y + float64(i)*cellSize + cellSize/2
68 px := min.X + cellSize/2
69 for j := range row {
70 v, ok := tree.Value(px, py)
71 if ok {
72 row[j] = v
73 }
74 px += cellSize
75 }
76 }
77 }
78
79 for n := runtime.NumCPU(); n >= 1; n-- {
80 wg.Add(1)
81 go rasterRow()
82 }
83
84 for i := 0; i < ycells; i++ {
85 rows <- i
86 }
87 close(rows)
88
89 wg.Wait()
90 log.Printf("info: Rastering took %v\n", time.Since(start))
91
92 start = time.Now()
93
94 tracer := contourmap.FromFloat64s(xcells+2, ycells+2, raster)
95
96 areas := make([]wkb.MultiPolygonGeom, len(heights))
97
98 // TODO: Check if this correct!
99 reprojX := common.Linear(1, min.X, float64(xcells+1), max.X)
100 reprojY := common.Linear(1, min.Y, float64(ycells+1), max.Y)
101
102 cnts := make(chan int)
103
104 doContours := func() {
105 defer wg.Done()
106 for hIdx := range cnts {
107 c := tracer.Contours(heights[hIdx])
108
109 if len(c) == 0 {
110 continue
111 }
112
113 // We need to bring it back to the
114 // none raster coordinate system.
115 a := make(wkb.MultiPolygonGeom, len(c))
116
117 for i, pl := range c {
118 shell := make(wkb.LinearRingGeom, len(pl))
119 for j, pt := range pl {
120 shell[j] = wkb.PointGeom{
121 X: reprojX(pt.X),
122 Y: reprojY(pt.Y),
123 }
124 }
125 a[i] = wkb.PolygonGeom{shell}
126 }
127
128 areas[hIdx] = a
129 }
130 }
131
132 for n := runtime.NumCPU(); n >= 1; n-- {
133 wg.Add(1)
134 go doContours()
135 }
136
137 for i := range heights {
138 cnts <- i
139 }
140 close(cnts)
141
142 wg.Wait()
143 log.Printf("info: Tracing areas took %v\n", time.Since(start))
144
145 return areas
146 }