Mercurial > gemma
changeset 4531:c9b6be8d81c8 iso-areas
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 | b74d98fc82a5 |
children | 769f723c2581 |
files | cmd/isoareas/main.go |
diffstat | 1 files changed, 111 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cmd/isoareas/main.go Tue Sep 24 17:19:52 2019 +0200 @@ -0,0 +1,111 @@ +// 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) 2019 by via donau +// – Österreichische Wasserstraßen-Gesellschaft mbH +// Software engineering by Intevation GmbH +// +// Author(s): +// * Sascha L. Teichmann <sascha.teichmann@intevation.de> + +package main + +import ( + "bufio" + "io" + "log" + "os" + "strconv" + "strings" + + "gemma.intevation.de/gemma/pkg/octree" +) + +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` + +func loadXYZ(r io.Reader) (octree.MultiPointZ, error) { + + scanner := bufio.NewScanner(r) + + points := make(octree.MultiPointZ, 0, 2000000) + + var x, y, z float64 + var err error + + for scanner.Scan() { + line := strings.TrimSpace(scanner.Text()) + if len(line) == 0 || strings.HasPrefix(line, "#") { + continue + } + parts := strings.SplitN(line, " ", 3) + if len(parts) != 3 { + continue + } + + if x, err = strconv.ParseFloat(parts[0], 64); err != nil { + return nil, err + } + if y, err = strconv.ParseFloat(parts[1], 64); err != nil { + return nil, err + } + if z, err = strconv.ParseFloat(parts[2], 64); err != nil { + return nil, err + } + points = append(points, octree.Vertex{X: x, Y: y, Z: z}) + } + + return points, nil +} + +func minMax(points []octree.Vertex) (octree.Vertex, octree.Vertex) { + if len(points) == 0 { + return octree.Vertex{}, octree.Vertex{} + } + + min, max := points[0], points[0] + + for _, v := range points { + min.Minimize(v) + max.Maximize(v) + } + + return min, max +} + +func check(err error) { + if err != nil { + log.Fatalf("error: %v\n", err) + } +} + +func main() { + + heights, err := octree.ParseClassBreaks(classBreaks) + check(err) + + xyz, err := loadXYZ(os.Stdin) + check(err) + + log.Printf("num vertices: %d\n", len(xyz)) + + min, max := minMax(xyz) + + log.Printf("(%.2f, %.2f, %.2f) - (%.2f, %.2f, %.2f)\n", + min.X, min.Y, min.Z, + max.X, max.Y, max.Z) + + heights = octree.ExtrapolateClassBreaks(heights, min.Z, max.Z) + + log.Println("classes:") + for i, v := range heights { + log.Printf("\t%d: %.3f\n", i+1, v) + } + + tri, err := octree.Triangulate(xyz) + check(err) + + _ = tri +}