Mercurial > gemma
view pkg/mesh/raster.go @ 5591:0011f50cf216 surveysperbottleneckid
Removed no longer used alternative api for surveys/ endpoint.
As bottlenecks in the summary for SR imports are now identified by
their id and no longer by the (not guarantied to be unique!) name,
there is no longer the need to request survey data by the name+date
tuple (which isn't reliable anyway). So the workaround was now
reversed.
author | Sascha Wilde <wilde@sha-bang.de> |
---|---|
date | Wed, 06 Apr 2022 13:30:29 +0200 |
parents | 5f47eeea988d |
children | 1222b777f51f |
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" ) type Raster struct { BBox Box2D CellSize float64 XCells int YCells int Cells []float64 } const noData = -math.MaxFloat64 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, } } 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() } 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() } 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 } 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 }