Mercurial > gemma
view cmd/isoareas/main.go @ 4556:04eba9dc917d iso-areas
Use colors from configuration instead of random rgb values.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Tue, 01 Oct 2019 17:37:54 +0200 |
parents | 1c5c6ffab886 |
children | 17cba8b447a6 |
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" "flag" "fmt" "io" "log" "math" "math/rand" "os" "runtime" "strconv" "strings" "sync" "time" svg "github.com/ajstarks/svgo" "github.com/fogleman/contourmap" "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 linear(x1, y1, x2, y2 float64) func(float64) float64 { // f(x1) = y1 // f(x2) = y2 // y1 = x1*a + b <=> b = y1 - x1*a // y2 = x2*a + b // y1 - y2 = a*(x1 - x2) // a = (y1-y2)/(x1 - x2) for x1 != x2 if x1 == x2 { return func(float64) float64 { return 0.5 * (y1 + y2) } } a := (y1 - y2) / (x1 - x2) b := y1 - x1*a return func(x float64) float64 { return x*a + b } } 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 := linear(minX, 0, maxX, width) py := 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 dumpPolygonsSVG( w io.Writer, min, max octree.Vertex, classPolygons [][]octree.LineStringZ, ) (err error) { ratio := (max.X - min.X) / (max.Y - min.Y) log.Printf("ratio: %.2f\n", ratio) const width = 50000 height := int(math.Ceil(width * ratio)) px := linear(min.X, 0, max.X, width) py := linear(min.Y, float64(height), max.Y, 0) out := bufio.NewWriter(w) defer func() { err = out.Flush() }() canvas := svg.New(out) canvas.Start(width, height) rnd := rand.New(rand.NewSource(42)) for _, polygons := range classPolygons { r := byte(rnd.Intn(256)) g := byte(rnd.Intn(256)) b := byte(rnd.Intn(256)) style := fmt.Sprintf("fill:#%02x%02x%02x", r, g, b) for _, polygon := range polygons { x := make([]int, len(polygon)) y := make([]int, len(polygon)) for i, p := range polygon { x[i] = int(math.Round(px(p.X))) y[i] = int(math.Round(py(p.Y))) } canvas.Polygon(x, y, style) } } canvas.End() return nil } func dumpArcsSVG( w io.Writer, min, max octree.Vertex, cuts [][]indexedArc, arcs []octree.LineStringZ, ) (err error) { ratio := (max.X - min.X) / (max.Y - min.Y) log.Printf("ratio: %.2f\n", ratio) const width = 50000 height := int(math.Ceil(width * ratio)) px := linear(min.X, 0, max.X, width) py := linear(min.Y, float64(height), max.Y, 0) out := bufio.NewWriter(w) defer func() { err = out.Flush() }() canvas := svg.New(out) canvas.Start(width, height) rnd := rand.New(rand.NewSource(42)) for _, cut := range cuts { usedArcs := map[int32]struct{}{} r := byte(rnd.Intn(256)) g := byte(rnd.Intn(256)) b := byte(rnd.Intn(256)) style := fmt.Sprintf("fill:none;stroke:#%02x%02x%02x", r, g, b) for i := range cut { idx := cut[i].arc if _, already := usedArcs[idx]; already { continue } usedArcs[idx] = struct{}{} arc := arcs[idx] x := make([]int, len(arc)) y := make([]int, len(arc)) for i, p := range arc { x[i] = int(math.Round(px(p.X))) y[i] = int(math.Round(py(p.Y))) } canvas.Polyline(x, y, style) } } canvas.End() return nil } 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)) */ }