Mercurial > gemma

comparison cmd/isoareas/main.go @ 4531:c9b6be8d81c8 iso-areas

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Started a branch to calculate the iso areas for a given triangulation.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Tue, 24 Sep 2019 17:19:52 +0200
parents
children 769f723c2581
comparison
equal deleted inserted replaced
4475:b74d98fc82a5 4531:c9b6be8d81c8
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) 2019 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 main
15
16 import (
17 "bufio"
18 "io"
19 "log"
20 "os"
21 "strconv"
22 "strings"
23
24 "gemma.intevation.de/gemma/pkg/octree"
25 )
26
27 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`
28
29 func loadXYZ(r io.Reader) (octree.MultiPointZ, error) {
30
31 scanner := bufio.NewScanner(r)
32
33 points := make(octree.MultiPointZ, 0, 2000000)
34
35 var x, y, z float64
36 var err error
37
38 for scanner.Scan() {
39 line := strings.TrimSpace(scanner.Text())
40 if len(line) == 0 || strings.HasPrefix(line, "#") {
41 continue
42 }
43 parts := strings.SplitN(line, " ", 3)
44 if len(parts) != 3 {
45 continue
46 }
47
48 if x, err = strconv.ParseFloat(parts[0], 64); err != nil {
49 return nil, err
50 }
51 if y, err = strconv.ParseFloat(parts[1], 64); err != nil {
52 return nil, err
53 }
54 if z, err = strconv.ParseFloat(parts[2], 64); err != nil {
55 return nil, err
56 }
57 points = append(points, octree.Vertex{X: x, Y: y, Z: z})
58 }
59
60 return points, nil
61 }
62
63 func minMax(points []octree.Vertex) (octree.Vertex, octree.Vertex) {
64 if len(points) == 0 {
65 return octree.Vertex{}, octree.Vertex{}
66 }
67
68 min, max := points[0], points[0]
69
70 for _, v := range points {
71 min.Minimize(v)
72 max.Maximize(v)
73 }
74
75 return min, max
76 }
77
78 func check(err error) {
79 if err != nil {
80 log.Fatalf("error: %v\n", err)
81 }
82 }
83
84 func main() {
85
86 heights, err := octree.ParseClassBreaks(classBreaks)
87 check(err)
88
89 xyz, err := loadXYZ(os.Stdin)
90 check(err)
91
92 log.Printf("num vertices: %d\n", len(xyz))
93
94 min, max := minMax(xyz)
95
96 log.Printf("(%.2f, %.2f, %.2f) - (%.2f, %.2f, %.2f)\n",
97 min.X, min.Y, min.Z,
98 max.X, max.Y, max.Z)
99
100 heights = octree.ExtrapolateClassBreaks(heights, min.Z, max.Z)
101
102 log.Println("classes:")
103 for i, v := range heights {
104 log.Printf("\t%d: %.3f\n", i+1, v)
105 }
106
107 tri, err := octree.Triangulate(xyz)
108 check(err)
109
110 _ = tri
111 }