Mercurial > gemma
view pkg/mesh/raster.go @ 5672:b1a10654bf0f
feat: Add Scale to map
author | Thomas Junk <thomas.junk@intevation.de> |
---|---|
date | Wed, 14 Feb 2024 16:56:28 +0100 |
parents | 1222b777f51f |
children | d2ccf6bb6940 |
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 mesh import ( "math" "runtime" "sync" "time" "github.com/fogleman/contourmap" "gemma.intevation.de/gemma/pkg/common" "gemma.intevation.de/gemma/pkg/log" "gemma.intevation.de/gemma/pkg/wkb" ) // Raster represents a 2D cell raster with a given extend. // The cells are equally sized in X and Y direction. type Raster struct { BBox Box2D CellSize float64 XCells int YCells int Cells []float64 } const noData = -math.MaxFloat64 // NewRaster creates a new Raster with the given extend // and cell size. func NewRaster(bbox Box2D, cellSize float64) *Raster { width, height := bbox.Size() log.Infof("raster extent: %.2f / %.2f", width, height) xCells := int(math.Ceil(width / cellSize)) yCells := int(math.Ceil(height / cellSize)) log.Infof("raster size: %d / %d\n", xCells, yCells) size := (xCells + 2) * (yCells + 2) cells := make([]float64, size) for i := range cells { cells[i] = noData } return &Raster{ BBox: bbox, CellSize: cellSize, XCells: xCells, YCells: yCells, Cells: cells, } } // Rasterize fills the raster with the evaluation data // received via the given eval function. The cell values // 4x oversampled, func (r *Raster) Rasterize(eval func(float64, float64) (float64, bool)) { var wg sync.WaitGroup rows := make(chan int) rasterRow := func() { defer wg.Done() quat := 0.25 * r.CellSize for i := range rows { pos := (i+1)*(r.XCells+2) + 1 row := r.Cells[pos : pos+r.XCells] py := r.BBox.Y1 + float64(i)*r.CellSize + r.CellSize/2 px := r.BBox.X1 + r.CellSize/2 y1 := py - quat y2 := py + quat for j := range row { var n int var sum float64 if v, ok := eval(px-quat, y1); ok { sum = v n = 1 } if v, ok := eval(px-quat, y2); ok { sum += v n++ } if v, ok := eval(px+quat, y1); ok { sum += v n++ } if v, ok := eval(px+quat, y2); ok { sum += v n++ } if n > 0 { row[j] = sum / float64(n) } px += r.CellSize } } } for n := runtime.NumCPU(); n >= 1; n-- { wg.Add(1) go rasterRow() } for i := 0; i < r.YCells; i++ { rows <- i } close(rows) wg.Wait() } // Diff updates the cell values of the raster with difference // of this raster and the values from the gien eval function. // The data will 4x over-sampled. func (r *Raster) Diff(eval func(float64, float64) (float64, bool)) { var wg sync.WaitGroup rows := make(chan int) rasterRow := func() { defer wg.Done() quat := 0.25 * r.CellSize for i := range rows { pos := (i+1)*(r.XCells+2) + 1 row := r.Cells[pos : pos+r.XCells] py := r.BBox.Y1 + float64(i)*r.CellSize + r.CellSize/2 px := r.BBox.X1 + r.CellSize/2 y1 := py - quat y2 := py + quat for j, old := range row { // only diff where values if old == noData { px += r.CellSize continue } var n int var sum float64 if v, ok := eval(px-quat, y1); ok { sum = v n = 1 } if v, ok := eval(px-quat, y2); ok { sum += v n++ } if v, ok := eval(px+quat, y1); ok { sum += v n++ } if v, ok := eval(px+quat, y2); ok { sum += v n++ } if n > 0 { row[j] -= sum / float64(n) } else { row[j] = noData } px += r.CellSize } } } for n := runtime.NumCPU(); n >= 1; n-- { wg.Add(1) go rasterRow() } for i := 0; i < r.YCells; i++ { rows <- i } close(rows) wg.Wait() } // ZExtent returns the z-range of the raster. // The last return value is false if the raster // only consists of no data values. func (r *Raster) ZExtent() (float64, float64, bool) { min, max := math.MaxFloat64, -math.MaxFloat64 for _, v := range r.Cells { if v == noData { continue } if v < min { min = v } if v > max { max = v } } return min, max, min != math.MaxFloat64 } // Trace generates contour lines of this raster given // the given class breaks. func (r *Raster) Trace(heights ClassBreaks) []wkb.MultiPolygonGeom { start := time.Now() tracer := contourmap.FromFloat64s(r.XCells+2, r.YCells+2, r.Cells) areas := make([]wkb.MultiPolygonGeom, len(heights)) reprojX := common.Linear(0.5, r.BBox.X1, 1.5, r.BBox.X1+r.CellSize) reprojY := common.Linear(0.5, r.BBox.Y1, 1.5, r.BBox.Y1+r.CellSize) var wg sync.WaitGroup cnts := make(chan int) doContours := func() { defer wg.Done() for hIdx := range cnts { if c := tracer.Contours(heights[hIdx]); len(c) > 0 { areas[hIdx] = buildMultipolygon(c, reprojX, reprojY) } } } for n := runtime.NumCPU(); n >= 1; n-- { wg.Add(1) go doContours() } for i := range heights { cnts <- i } close(cnts) wg.Wait() log.Infof("Tracing areas took %v\n", time.Since(start)) return areas } type contour []contourmap.Point type bboxNode struct { box Box2D cnt contour children []*bboxNode } func (cnt contour) contains(o contour) bool { return contains(cnt, o[0].X, o[0].Y) || contains(cnt, o[len(o)/2].X, o[len(o)/2].Y) } func (cnt contour) length() int { return len(cnt) } func (cnt contour) point(i int) (float64, float64) { return cnt[i].X, cnt[i].Y } func (cnt contour) bbox() Box2D { minX, minY := math.MaxFloat64, math.MaxFloat64 maxX, maxY := -math.MaxFloat64, -math.MaxFloat64 for _, p := range cnt { if p.X < minX { minX = p.X } if p.X > maxX { maxX = p.X } if p.Y < minY { minY = p.Y } if p.Y > maxY { maxY = p.Y } } return Box2D{X1: minX, X2: maxX, Y1: minY, Y2: maxY} } func (bn *bboxNode) insert(cnt contour, box Box2D) { // check if children are inside new var nr *bboxNode for i, r := range bn.children { if r.box.Inside(box) && cnt.contains(r.cnt) { if nr == nil { nr = &bboxNode{box: box, cnt: cnt} } nr.children = append(nr.children, r) bn.children[i] = nil } } // we have a new child if nr != nil { // compact the list for i := len(bn.children) - 1; i >= 0; i-- { if bn.children[i] == nil { if i < len(bn.children)-1 { copy(bn.children[i:], bn.children[i+1:]) } bn.children[len(bn.children)-1] = nil bn.children = bn.children[:len(bn.children)-1] } } bn.children = append(bn.children, nr) return } // check if new is inside an old for _, r := range bn.children { if box.Inside(r.box) && r.cnt.contains(cnt) { r.insert(cnt, box) return } } // its a new child node. nr = &bboxNode{box: box, cnt: cnt} bn.children = append(bn.children, nr) } func (bn *bboxNode) insertRoot(cnt contour) { bn.insert(cnt, cnt.bbox()) } type bboxOutFunc func(contour, []contour) func (bn *bboxNode) generate(out bboxOutFunc) { var grands []*bboxNode holes := make([]contour, len(bn.children)) for i, ch := range bn.children { holes[i] = ch.cnt grands = append(grands, ch.children...) } out(bn.cnt, holes) // the grand children are new polygons. for _, grand := range grands { grand.generate(out) } } func (bn *bboxNode) generateRoot(out bboxOutFunc) { for _, r := range bn.children { r.generate(out) } } func buildMultipolygon( cnts []contourmap.Contour, reprojX, reprojY func(float64) float64, ) wkb.MultiPolygonGeom { var forest bboxNode for _, cnt := range cnts { forest.insertRoot(contour(cnt)) } //log.Debugf("cnts: %d roots: %d\n", len(cnts), len(bf.roots)) var mp wkb.MultiPolygonGeom out := func(sh contour, hls []contour) { polygon := make(wkb.PolygonGeom, 1+len(hls)) // Handle shell shell := make(wkb.LinearRingGeom, len(sh)) for j, pt := range sh { shell[j] = wkb.PointGeom{ X: reprojX(pt.X), Y: reprojY(pt.Y), } } if shell.CCW() { shell.Reverse() } polygon[0] = shell // handle holes for i, hl := range hls { hole := make(wkb.LinearRingGeom, len(hl)) for j, pt := range hl { hole[j] = wkb.PointGeom{ X: reprojX(pt.X), Y: reprojY(pt.Y), } } if !hole.CCW() { hole.Reverse() } polygon[1+i] = hole } mp = append(mp, polygon) } forest.generateRoot(out) return mp }