Mercurial > gemma
comparison 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 |
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 } |