Mercurial > gemma
view cmd/isoareas/main.go @ 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 | |
children | 769f723c2581 |
line wrap: on
line source
// 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 }