Mercurial > gemma

comparison pkg/imports/sr.go @ 3600:6590733e2ebc

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
SR import: If in single beam mode and no boundary is given eliminate triangles which have at least an edge with is longer than 98% of the overall edge lengths.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Tue, 04 Jun 2019 18:37:15 +0200
parents eeeb7bf14217
children 6248a4bc10c7
comparison
equal deleted inserted replaced
3599:73add36de413 3600:6590733e2ebc
26 "io" 26 "io"
27 "math" 27 "math"
28 "os" 28 "os"
29 "path" 29 "path"
30 "path/filepath" 30 "path/filepath"
31 "sort"
31 "strconv" 32 "strconv"
32 "strings" 33 "strings"
33 "time" 34 "time"
34 35
35 shp "github.com/jonas-p/go-shp" 36 shp "github.com/jonas-p/go-shp"
37 "gonum.org/v1/gonum/stat"
36 38
37 "gemma.intevation.de/gemma/pkg/common" 39 "gemma.intevation.de/gemma/pkg/common"
38 "gemma.intevation.de/gemma/pkg/misc" 40 "gemma.intevation.de/gemma/pkg/misc"
39 "gemma.intevation.de/gemma/pkg/models" 41 "gemma.intevation.de/gemma/pkg/models"
40 "gemma.intevation.de/gemma/pkg/octree" 42 "gemma.intevation.de/gemma/pkg/octree"
174 AND grwl.validity = bns.gauge_validity 176 AND grwl.validity = bns.gauge_validity
175 WHERE bns.objnam = $1 AND grwl.depth_reference like 'LDC%' 177 WHERE bns.objnam = $1 AND grwl.depth_reference like 'LDC%'
176 ` 178 `
177 179
178 reprojectPointsSingleBeamSQL = ` 180 reprojectPointsSingleBeamSQL = `
179 SELECT ST_AsBinary( 181 SELECT
180 ST_Transform( 182 ST_AsBinary(
181 ST_GeomFromWKB($1, $2::integer), 183 ST_Transform(
182 best_utm(CAST(ST_GeomFromWKB($1, $2::integer) AS geography))), 184 ST_GeomFromWKB($1, $2::integer),
183 best_utm(CAST(ST_GeomFromWKB($1, $2::integer) AS geography)) 185 best_utm(CAST(ST_Transform(ST_GeomFromWKB($1, $2::integer), 4326) AS geography)))),
186 best_utm(CAST(ST_Transform(ST_GeomFromWKB($1, $2::integer), 4326) AS geography))
184 ` 187 `
185 ) 188 )
186 189
187 func (sr *SoundingResult) isSingleBeam() bool { 190 func (sr *SoundingResult) isSingleBeam() bool {
188 return sr.SingleBeam != nil && *sr.SingleBeam 191 return sr.SingleBeam != nil && *sr.SingleBeam
203 return nil, err 206 return nil, err
204 } 207 }
205 defer z.Close() 208 defer z.Close()
206 209
207 feedback.Info("Looking for 'meta.json'") 210 feedback.Info("Looking for 'meta.json'")
211
208 mf := common.FindInZIP(z, "meta.json") 212 mf := common.FindInZIP(z, "meta.json")
209 if mf == nil && !sr.completeOverride() { 213 if mf == nil && !sr.completeOverride() {
210 return nil, errors.New("Cannot find 'meta.json'") 214 return nil, errors.New("Cannot find 'meta.json'")
211 } 215 }
212 216
340 m.EPSG, 344 m.EPSG,
341 ).Scan(&reproj, &epsg); err != nil { 345 ).Scan(&reproj, &epsg); err != nil {
342 return nil, err 346 return nil, err
343 } 347 }
344 348
349 if err := xyz.FromWKB(reproj); err != nil {
350 return nil, err
351 }
352
345 feedback.Info("Reprojecting points to EPSG %d took %v.", 353 feedback.Info("Reprojecting points to EPSG %d took %v.",
346 epsg, time.Since(start)) 354 epsg, time.Since(start))
355 feedback.Info("Number of reprojected points: %d", len(xyz))
356
357 start = time.Now()
358 tri, err := octree.Triangulate(xyz)
359 if err != nil {
360 return nil, err
361 }
362 feedback.Info("Triangulation took %v.", time.Since(start))
363 feedback.Info("Number triangles: %d.", len(tri.Triangles)/3)
364
365 if boundary == nil {
366 feedback.Info("No boundary given.")
367 feedback.Info("Eliminate triangles with long edges.")
368
369 lens := make([]float64, len(tri.Triangles)/3)
370 for i := 0; i < len(tri.Triangles); i += 3 {
371 vs := tri.Triangles[i : i+3]
372 var max float64
373 for j := 0; j < 3; j++ {
374 k := (j + 1) % 3
375 if l := tri.Points[vs[j]].SquaredDistance2D(tri.Points[vs[k]]); l > max {
376 max = l
377 }
378 }
379 lens[i/3] = max
380 }
381 qs := make([]float64, len(lens))
382 copy(qs, lens)
383 sort.Float64s(qs)
384 q := stat.Quantile(0.98, stat.LinInterp, qs, nil)
385 feedback.Info("98%% quantile of the edge lengths: %.2f\n", math.Sqrt(q))
386
387 ts := tri.Triangles
388
389 for i := len(lens) - 1; i >= 0; i-- {
390 if lens[i] >= q {
391 if i != len(lens)-1 {
392 copy(ts[i*3:], ts[(i+1)*3:])
393 }
394 ts = ts[:len(ts)-3]
395 }
396 }
397 feedback.Info("Remaining triangles: %d (%.2f%%)\n",
398 len(ts)/3, float64(len(ts))/float64(len(tri.Triangles))*100.0)
399 }
347 400
348 // TODO: Implement me! 401 // TODO: Implement me!
349 return nil, errors.New("Not implemented, yet!") 402 return nil, errors.New("Not implemented, yet!")
350 } 403 }
351 404