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
}