Mercurial > gemma

comparison cmd/isoareas/main.go @ 4552:7ed5a4a94efc iso-areas

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Used fogleman's contourmap as a tracing alternative.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Tue, 01 Oct 2019 13:18:16 +0200
parents b5b23b6d76f1
children 4c476d65d1bb
comparison
equal deleted inserted replaced
4551:b5b23b6d76f1 4552:7ed5a4a94efc
13 13
14 package main 14 package main
15 15
16 import ( 16 import (
17 "bufio" 17 "bufio"
18 "flag"
18 "fmt" 19 "fmt"
19 "io" 20 "io"
20 "log" 21 "log"
21 "math" 22 "math"
22 "math/rand" 23 "math/rand"
23 "os" 24 "os"
25 "runtime"
24 "strconv" 26 "strconv"
25 "strings" 27 "strings"
28 "sync"
26 "time" 29 "time"
27 30
28 svg "github.com/ajstarks/svgo" 31 svg "github.com/ajstarks/svgo"
32 "github.com/fogleman/contourmap"
29 33
30 "gemma.intevation.de/gemma/pkg/octree" 34 "gemma.intevation.de/gemma/pkg/octree"
31 ) 35 )
32 36
33 const classBreaks = `1:#ff00dd,1.5,1.7,1.9,2.1,2.3,2.5:#f25f20,2.7,2.9,3.1:#f7e40e,3.3,3.5,4:#8ad51a,4.5,5,5.5,6,6.5,7:#1414ff` 37 const classBreaks = `1:#ff00dd,1.5,1.7,1.9,2.1,2.3,2.5:#f25f20,2.7,2.9,3.1:#f7e40e,3.3,3.5,4:#8ad51a,4.5,5,5.5,6,6.5,7:#1414ff`
62 } 66 }
63 points = append(points, octree.Vertex{X: x, Y: y, Z: z}) 67 points = append(points, octree.Vertex{X: x, Y: y, Z: z})
64 } 68 }
65 69
66 return points, nil 70 return points, nil
67 }
68
69 func minMax(points []octree.Vertex) (octree.Vertex, octree.Vertex) {
70 if len(points) == 0 {
71 return octree.Vertex{}, octree.Vertex{}
72 }
73
74 min, max := points[0], points[0]
75
76 for _, v := range points {
77 min.Minimize(v)
78 max.Maximize(v)
79 }
80
81 return min, max
82 } 71 }
83 72
84 func check(err error) { 73 func check(err error) {
85 if err != nil { 74 if err != nil {
86 log.Fatalf("error: %v\n", err) 75 log.Fatalf("error: %v\n", err)
104 a := (y1 - y2) / (x1 - x2) 93 a := (y1 - y2) / (x1 - x2)
105 b := y1 - x1*a 94 b := y1 - x1*a
106 return func(x float64) float64 { 95 return func(x float64) float64 {
107 return x*a + b 96 return x*a + b
108 } 97 }
98 }
99
100 func dumpContoursSVG(
101 w io.Writer,
102 contours [][]contourmap.Contour,
103 ) (err error) {
104
105 var (
106 minX = math.MaxFloat64
107 minY = math.MaxFloat64
108 maxX = -math.MaxFloat64
109 maxY = -math.MaxFloat64
110 )
111
112 for _, c := range contours {
113 for _, r := range c {
114 for _, p := range r {
115 if p.X < minX {
116 minX = p.X
117 }
118 if p.Y < minY {
119 minY = p.Y
120 }
121 if p.X > maxX {
122 maxX = p.X
123 }
124 if p.Y > maxY {
125 maxY = p.Y
126 }
127 }
128 }
129 }
130 ratio := (maxX - minX) / (maxY - minY)
131
132 log.Printf("ratio: %.2f\n", ratio)
133
134 const width = 50000
135 height := int(math.Ceil(width * ratio))
136
137 px := linear(minX, 0, maxX, width)
138 py := linear(minY, float64(height), maxY, 0)
139
140 out := bufio.NewWriter(w)
141 defer func() { err = out.Flush() }()
142
143 canvas := svg.New(out)
144 canvas.Start(width, height)
145
146 rnd := rand.New(rand.NewSource(42))
147
148 for _, c := range contours {
149
150 r := byte(rnd.Intn(256))
151 g := byte(rnd.Intn(256))
152 b := byte(rnd.Intn(256))
153 style := fmt.Sprintf("fill:#%02x%02x%02x", r, g, b)
154
155 for _, r := range c {
156 x := make([]int, len(r))
157 y := make([]int, len(r))
158 for i, p := range r {
159 x[i] = int(math.Round(px(p.X)))
160 y[i] = int(math.Round(py(p.Y)))
161 }
162
163 canvas.Polygon(x, y, style)
164 }
165 }
166
167 canvas.End()
168 return
109 } 169 }
110 170
111 func dumpPolygonsSVG( 171 func dumpPolygonsSVG(
112 w io.Writer, 172 w io.Writer,
113 min, max octree.Vertex, 173 min, max octree.Vertex,
216 return nil 276 return nil
217 } 277 }
218 278
219 func main() { 279 func main() {
220 280
281 cellSize := float64(0.5)
282
283 flag.Float64Var(&cellSize, "cellsize", 0.5, "width/height of raster cell [m]")
284
285 flag.Parse()
286
221 heights, err := octree.ParseClassBreaks(classBreaks) 287 heights, err := octree.ParseClassBreaks(classBreaks)
222 check(err) 288 check(err)
223 log.Printf("num classes: %d\n", len(heights)) 289 log.Printf("num classes: %d\n", len(heights))
224 290
225 start := time.Now() 291 start := time.Now()
227 check(err) 293 check(err)
228 log.Printf("loading took %v\n", time.Since(start)) 294 log.Printf("loading took %v\n", time.Since(start))
229 295
230 log.Printf("num vertices: %d\n", len(xyz)) 296 log.Printf("num vertices: %d\n", len(xyz))
231 297
232 min, max := minMax(xyz) 298 start = time.Now()
299 tri, err := octree.Triangulate(xyz)
300 check(err)
301 log.Printf("triangulation took %v\n", time.Since(start))
302 tooLongEdge := tri.EstimateTooLong()
303 log.Printf("Eliminate triangles with edges longer than %.2f meters.", tooLongEdge)
304
305 start = time.Now()
306 polygon, removed := tri.ConcaveHull(tooLongEdge)
307
308 polygonArea := polygon.Area()
309 if polygonArea < 0.0 { // counter clockwise
310 polygonArea = -polygonArea
311 polygon.Reverse()
312 }
313
314 var clippingPolygon octree.Polygon
315 clippingPolygon.FromLineStringZ(polygon)
316 clippingPolygon.Indexify()
317
318 tin := tri.Tin()
319 // tin.EPSG = epsg
320
321 var str octree.STRTree
322 str.Build(tin)
323
324 removed = str.Clip(&clippingPolygon)
325
326 builder := octree.NewBuilder(tin)
327 builder.Build(removed)
328
329 tree := builder.Tree()
330
331 min, max := tin.Min, tin.Max
233 332
234 log.Printf("(%.2f, %.2f, %.2f) - (%.2f, %.2f, %.2f)\n", 333 log.Printf("(%.2f, %.2f, %.2f) - (%.2f, %.2f, %.2f)\n",
235 min.X, min.Y, min.Z, 334 min.X, min.Y, min.Z,
236 max.X, max.Y, max.Z) 335 max.X, max.Y, max.Z)
237 336
238 start = time.Now() 337 width := max.X - min.X
239 tri, err := octree.Triangulate(xyz) 338 height := max.Y - min.Y
240 check(err) 339
241 log.Printf("triangulation took %v\n", time.Since(start)) 340 log.Printf("width/height: %.2f / %.2f\n", width, height)
242 341
243 log.Printf("number of triangles: %d\n", len(tri.Triangles)/3) 342 xcells := int(math.Ceil(width / cellSize))
244 343 ycells := int(math.Ceil(height / cellSize))
245 start = time.Now() 344
246 polygons := intersectTriangles(tri, heights, min, max) 345 log.Printf("raster size: (%d, %d)\n", xcells, ycells)
247 log.Printf("intersecting triangles took %v.\n", time.Since(start)) 346
248 347 raster := make([]float64, xcells*ycells)
249 check(dumpPolygonsSVG(os.Stdout, min, max, polygons)) 348
250 } 349 const closed = -math.MaxFloat64
350 for i := range raster {
351 raster[i] = closed
352 }
353
354 start = time.Now()
355
356 var wg sync.WaitGroup
357
358 rows := make(chan int)
359
360 rasterRow := func() {
361 defer wg.Done()
362 for i := range rows {
363 pos := i * xcells
364 row := raster[pos : pos+xcells]
365 //log.Printf("len: %d\n", len(row))
366 py := min.Y + float64(i)*cellSize + cellSize/2
367 px := min.X + cellSize/2
368 for j := range row {
369 v, ok := tree.Value(px, py)
370 if ok {
371 row[j] = v
372 }
373 px += cellSize
374 }
375 }
376 }
377
378 start = time.Now()
379
380 for n := runtime.NumCPU(); n >= 1; n-- {
381 wg.Add(1)
382 go rasterRow()
383 }
384
385 for i := 0; i < ycells; i++ {
386 rows <- i
387 }
388 close(rows)
389
390 wg.Wait()
391
392 log.Printf("Rasterizing took %v.\n", time.Since(start))
393
394 start = time.Now()
395 cm := contourmap.FromFloat64s(xcells, ycells, raster).Closed()
396
397 start = time.Now()
398
399 contours := make([][]contourmap.Contour, len(heights))
400
401 cnts := make(chan int)
402
403 doContours := func() {
404 defer wg.Done()
405 for i := range cnts {
406 contours[i] = cm.Contours(heights[i])
407 }
408 }
409
410 for n := runtime.NumCPU(); n >= 1; n-- {
411 wg.Add(1)
412 go doContours()
413 }
414
415 for i := range heights {
416 cnts <- i
417 }
418 close(cnts)
419
420 wg.Wait()
421
422 log.Printf("Calculating contours took %v.\n", time.Since(start))
423 check(dumpContoursSVG(os.Stdout, contours))
424
425 /*
426
427 start = time.Now()
428 polygons := intersectTriangles(tri, heights, min, max)
429 log.Printf("intersecting triangles took %v.\n", time.Since(start))
430
431 check(dumpPolygonsSVG(os.Stdout, min, max, polygons))
432 */
433 }