Mercurial > gemma
changeset 4588:7844a3a630eb iso-areas
Removed obsolete dev tool for experimenting with iso areas.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Tue, 08 Oct 2019 22:42:21 +0200 |
parents | ffa3148b95c5 |
children | acd802a76b93 |
files | cmd/isoareas/main.go go.mod go.sum |
diffstat | 3 files changed, 0 insertions(+), 313 deletions(-) [+] |
line wrap: on
line diff
--- a/cmd/isoareas/main.go Tue Oct 08 22:36:54 2019 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,310 +0,0 @@ -// 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" - "flag" - "fmt" - "io" - "log" - "math" - "os" - "runtime" - "strconv" - "strings" - "sync" - "time" - - svg "github.com/ajstarks/svgo" - "github.com/fogleman/contourmap" - - "gemma.intevation.de/gemma/pkg/common" - "gemma.intevation.de/gemma/pkg/models" - "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 check(err error) { - if err != nil { - log.Fatalf("error: %v\n", err) - } -} - -func dumpContoursSVG( - w io.Writer, - heights []float64, - colors models.ColorValues, - contours [][]contourmap.Contour, -) (err error) { - - var ( - minX = math.MaxFloat64 - minY = math.MaxFloat64 - maxX = -math.MaxFloat64 - maxY = -math.MaxFloat64 - ) - - for _, c := range contours { - for _, r := range c { - for _, p := range r { - if p.X < minX { - minX = p.X - } - if p.Y < minY { - minY = p.Y - } - if p.X > maxX { - maxX = p.X - } - if p.Y > maxY { - maxY = p.Y - } - } - } - } - ratio := (maxX - minX) / (maxY - minY) - - log.Printf("ratio: %.2f\n", ratio) - - const width = 50000 - height := int(math.Ceil(width * ratio)) - - px := common.Linear(minX, 0, maxX, width) - py := common.Linear(minY, float64(height), maxY, 0) - - out := bufio.NewWriter(w) - defer func() { err = out.Flush() }() - - canvas := svg.New(out) - canvas.Start(width, height) - - for j, c := range contours { - - h := heights[j] - col := colors.Clip(h) - - style := fmt.Sprintf("fill:#%02x%02x%02x", col.R, col.G, col.B) - - for _, r := range c { - x := make([]int, len(r)) - y := make([]int, len(r)) - for i, p := range r { - x[i] = int(math.Round(px(p.X))) - y[i] = int(math.Round(py(p.Y))) - } - - canvas.Polygon(x, y, style) - } - } - - canvas.End() - return -} - -func main() { - - cellSize := float64(0.5) - - flag.Float64Var(&cellSize, "cellsize", 0.5, "width/height of raster cell [m]") - - flag.Parse() - - colors, err := models.ParseColorValues(classBreaks) - check(err) - - heights := colors.Heights() - - log.Printf("num classes: %d\n", len(heights)) - start := time.Now() - xyz, err := loadXYZ(os.Stdin) - check(err) - log.Printf("loading took %v\n", time.Since(start)) - - log.Printf("num vertices: %d\n", len(xyz)) - - start = time.Now() - tri, err := octree.Triangulate(xyz) - check(err) - log.Printf("triangulation took %v\n", time.Since(start)) - tooLongEdge := tri.EstimateTooLong() - log.Printf("Eliminate triangles with edges longer than %.2f meters.", tooLongEdge) - - start = time.Now() - polygon, removed := tri.ConcaveHull(tooLongEdge) - - polygonArea := polygon.Area() - if polygonArea < 0.0 { // counter clockwise - polygonArea = -polygonArea - polygon.Reverse() - } - - var clippingPolygon octree.Polygon - clippingPolygon.FromLineStringZ(polygon) - clippingPolygon.Indexify() - - tin := tri.Tin() - // tin.EPSG = epsg - - var str octree.STRTree - str.Build(tin) - - removed = str.Clip(&clippingPolygon) - - builder := octree.NewBuilder(tin) - builder.Build(removed) - - tree := builder.Tree() - - min, max := tin.Min, tin.Max - - 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) - - width := max.X - min.X - height := max.Y - min.Y - - log.Printf("width/height: %.2f / %.2f\n", width, height) - - xcells := int(math.Ceil(width / cellSize)) - ycells := int(math.Ceil(height / cellSize)) - - log.Printf("raster size: (%d, %d)\n", xcells, ycells) - - // Add border for closing - raster := make([]float64, (xcells+2)*(ycells+2)) - - const closed = -math.MaxFloat64 - for i := range raster { - raster[i] = closed - } - - start = time.Now() - - var wg sync.WaitGroup - - rows := make(chan int) - - rasterRow := func() { - defer wg.Done() - for i := range rows { - pos := (i+1)*(xcells+2) + 1 - row := raster[pos : pos+xcells] - //log.Printf("len: %d\n", len(row)) - py := min.Y + float64(i)*cellSize + cellSize/2 - px := min.X + cellSize/2 - for j := range row { - v, ok := tree.Value(px, py) - if ok { - row[j] = v - } - px += cellSize - } - } - } - - start = time.Now() - - for n := runtime.NumCPU(); n >= 1; n-- { - wg.Add(1) - go rasterRow() - } - - for i := 0; i < ycells; i++ { - rows <- i - } - close(rows) - - wg.Wait() - - log.Printf("Rasterizing took %v.\n", time.Since(start)) - - start = time.Now() - cm := contourmap.FromFloat64s(xcells+2, ycells+2, raster) - - start = time.Now() - - contours := make([][]contourmap.Contour, len(heights)) - - cnts := make(chan int) - - doContours := func() { - defer wg.Done() - for i := range cnts { - contours[i] = cm.Contours(heights[i]) - } - } - - for n := runtime.NumCPU(); n >= 1; n-- { - wg.Add(1) - go doContours() - } - - for i := range heights { - cnts <- i - } - close(cnts) - - wg.Wait() - - log.Printf("Calculating contours took %v.\n", time.Since(start)) - check(dumpContoursSVG(os.Stdout, heights, colors, contours)) - - /* - - start = time.Now() - polygons := intersectTriangles(tri, heights, min, max) - log.Printf("intersecting triangles took %v.\n", time.Since(start)) - - check(dumpPolygonsSVG(os.Stdout, min, max, polygons)) - */ -}
--- a/go.mod Tue Oct 08 22:36:54 2019 +0200 +++ b/go.mod Tue Oct 08 22:42:21 2019 +0200 @@ -3,7 +3,6 @@ go 1.13 require ( - github.com/ajstarks/svgo v0.0.0-20180226025133-644b8db467af github.com/cockroachdb/apd v1.1.0 // indirect github.com/etcd-io/bbolt v1.3.3 github.com/fogleman/contourmap v0.0.0-20190814184649-9f61d36c4199
--- a/go.sum Tue Oct 08 22:36:54 2019 +0200 +++ b/go.sum Tue Oct 08 22:42:21 2019 +0200 @@ -64,8 +64,6 @@ github.com/jackc/fake v0.0.0-20150926172116-812a484cc733/go.mod h1:WrMFNQdiFJ80sQsxDoMokWK1W5TQtxBFNpzWTD84ibQ= github.com/jackc/pgx v3.6.0+incompatible h1:bJeo4JdVbDAW8KB2m8XkFeo8CPipREoG37BwEoKGz+Q= github.com/jackc/pgx v3.6.0+incompatible/go.mod h1:0ZGrqGqkRlliWnWB4zKnWtjbSWbGkVEFm4TeybAXq+I= -github.com/jonas-p/go-shp v0.1.1 h1:LY81nN67DBCz6VNFn2kS64CjmnDo9IP8rmSkTvhO9jE= -github.com/jonas-p/go-shp v0.1.1/go.mod h1:MRIhyxDQ6VVp0oYeD7yPGr5RSTNScUFKCDsI5DR7PtI= github.com/jonas-p/go-shp v0.1.2-0.20190401125246-9fd306ae10a6 h1:h5O7ee4tlSPVjdC75eSLX7jXZiHftthuHio/GtrhaSM= github.com/jonas-p/go-shp v0.1.2-0.20190401125246-9fd306ae10a6/go.mod h1:MRIhyxDQ6VVp0oYeD7yPGr5RSTNScUFKCDsI5DR7PtI= github.com/jonboulle/clockwork v0.1.0/go.mod h1:Ii8DK3G1RaLaWxj9trq07+26W01tbo22gdxWY5EU2bo=