Mercurial > gemma

comparison pkg/imports/sr.go @ 3655:a6c671abbc35 single-beam

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Started with cleaning up the single beam import.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Thu, 13 Jun 2019 18:33:54 +0200
parents 26e02fa5b992
children 1c3df921361d
comparison
equal deleted inserted replaced
3654:26e02fa5b992 3655:a6c671abbc35
22 "database/sql" 22 "database/sql"
23 "encoding/hex" 23 "encoding/hex"
24 "errors" 24 "errors"
25 "fmt" 25 "fmt"
26 "io" 26 "io"
27 "log"
28 "math" 27 "math"
29 "os" 28 "os"
30 "path" 29 "path"
31 "path/filepath" 30 "path/filepath"
32 "runtime"
33 "strconv" 31 "strconv"
34 "strings" 32 "strings"
35 "sync"
36 "time" 33 "time"
37 34
38 shp "github.com/jonas-p/go-shp" 35 shp "github.com/jonas-p/go-shp"
39 36
40 "gemma.intevation.de/gemma/pkg/common" 37 "gemma.intevation.de/gemma/pkg/common"
148 145
149 repairBoundarySQL = ` 146 repairBoundarySQL = `
150 SELECT 147 SELECT
151 ST_AsBinary(ST_Buffer(ST_MakeValid(ST_GeomFromWKB($1, $2::integer)), 0.0)), 148 ST_AsBinary(ST_Buffer(ST_MakeValid(ST_GeomFromWKB($1, $2::integer)), 0.0)),
152 ST_AsBinary(ST_Buffer(ST_MakeValid(ST_GeomFromWKB($1, $2::integer)), 0.1))` 149 ST_AsBinary(ST_Buffer(ST_MakeValid(ST_GeomFromWKB($1, $2::integer)), 0.1))`
153
154 whatsWrongSQL = `
155 SELECT
156 ST_IsValidReason(
157 CAST(ST_MakeValid(ST_Transform(ST_GeomFromWKB($1, $2::integer), 4326))::geography AS geometry)),
158 ST_IsValid(
159 CAST(ST_MakeValid(ST_Transform(ST_GeomFromWKB($1, $2::integer), 4326))::geography AS geometry))
160 `
161 150
162 insertContourSQL = ` 151 insertContourSQL = `
163 INSERT INTO waterway.sounding_results_contour_lines ( 152 INSERT INTO waterway.sounding_results_contour_lines (
164 sounding_result_id, 153 sounding_result_id,
165 height, 154 height,
382 feedback.Info("Number triangles: %d.", len(tri.Triangles)/3) 371 feedback.Info("Number triangles: %d.", len(tri.Triangles)/3)
383 372
384 var clippingPolygon octree.Polygon 373 var clippingPolygon octree.Polygon
385 374
386 if boundary == nil { 375 if boundary == nil {
387 feedback.Info("No boundary given.") 376 feedback.Info("No boundary given. Calulating from XYZ data.")
377 feedback.Info("Eliminate triangles with edges longer than %.2f meters.", tooLongEdge)
388 378
389 polygon, removed := tri.ConcaveHull(tooLongEdge) 379 polygon, removed := tri.ConcaveHull(tooLongEdge)
390 380
391 polygonArea := polygon.Area() 381 polygonArea := polygon.Area()
392 if polygonArea < 0.0 { // counter clockwise 382 if polygonArea < 0.0 { // counter clockwise
393 polygonArea = -polygonArea 383 polygonArea = -polygonArea
394 polygon.Reverse() 384 polygon.Reverse()
395 } 385 }
396 386
397 clippingPolygon.FromLineStringZ(polygon) 387 clippingPolygon.FromLineStringZ(polygon)
398
399 var wrong string
400 var valid bool
401 if err := tx.QueryRowContext(
402 ctx,
403 whatsWrongSQL,
404 clippingPolygon.AsWKB(),
405 epsg,
406 ).Scan(&wrong, &valid); err != nil {
407 return nil, err
408 }
409
410 log.Printf("!!!!!!!!: whats wrong (%t): %s\n", valid, wrong)
411 388
412 var repaired, buffered []byte 389 var repaired, buffered []byte
413 390
414 if err := tx.QueryRowContext( 391 if err := tx.QueryRowContext(
415 ctx, 392 ctx,
432 builder := octree.NewBuilder(firstTin) 409 builder := octree.NewBuilder(firstTin)
433 builder.Build(removed) 410 builder.Build(removed)
434 411
435 firstTree := builder.Tree() 412 firstTree := builder.Tree()
436 413
437 bbox := polygon.BBox() 414 feedback.Info("Boundary area: %.2fm²", polygonArea)
438 bboxArea := bbox.Area()
439
440 log.Printf("bbox area: %.2f\n", bboxArea)
441 log.Printf("polygon area: %.2f\n", polygonArea)
442 415
443 numPoints := int(math.Ceil(polygonArea * pointsPerSquareMeter)) 416 numPoints := int(math.Ceil(polygonArea * pointsPerSquareMeter))
444 417
445 log.Printf("random points for density of ~%d points/m²: %d\n", 418 feedback.Info("Generating %d random points for an average density of ~%d points/m².",
446 pointsPerSquareMeter, numPoints) 419 numPoints, pointsPerSquareMeter)
447 420
448 start := time.Now() 421 start := time.Now()
449 422
450 var wg sync.WaitGroup
451
452 jobs := make(chan int)
453 out := make(chan []octree.Vertex)
454 done := make(chan struct{})
455
456 cpus := runtime.NumCPU()
457
458 free := make(chan []octree.Vertex, cpus)
459
460 for i, n := 0, cpus; i < n; i++ {
461 wg.Add(1)
462 go func() {
463 defer wg.Done()
464 xRange := common.Random(bbox.X1, bbox.X2)
465 yRange := common.Random(bbox.Y1, bbox.Y2)
466
467 for size := range jobs {
468 var vertices []octree.Vertex
469 select {
470 case vertices = <-free:
471 default:
472 vertices = make([]octree.Vertex, 0, 1000)
473 }
474 for len(vertices) < size {
475 x, y := xRange(), yRange()
476 if z, ok := firstTree.Value(x, y); ok {
477 vertices = append(vertices, octree.Vertex{X: x, Y: y, Z: z})
478 }
479 }
480 out <- vertices
481 }
482 }()
483 }
484
485 generated := make(octree.LineStringZ, 0, numPoints+len(polygon)-1) 423 generated := make(octree.LineStringZ, 0, numPoints+len(polygon)-1)
486 424
487 go func() { 425 firstTree.GenerateRandomVertices(numPoints, func(vertices []octree.Vertex) {
488 defer close(done) 426 generated = append(generated, vertices...)
489 for vertices := range out { 427 })
490 generated = append(generated, vertices...) 428
491 select { 429 feedback.Info("Generating %d points took %v.", len(generated), time.Since(start))
492 case free <- vertices[:0]:
493 default:
494 }
495 }
496 }()
497
498 for remain := numPoints; remain > 0; {
499 if remain > 1000 {
500 jobs <- 1000
501 remain -= 1000
502 } else {
503 jobs <- remain
504 remain = 0
505 }
506 }
507 close(jobs)
508 wg.Wait()
509 close(out)
510 <-done
511
512 log.Printf("Generating %d points took %v\n", len(generated), time.Since(start))
513 430
514 // Add the boundary to new point cloud. 431 // Add the boundary to new point cloud.
515 generated = append(generated, polygon[:len(polygon)-1]...) 432 generated = append(generated, polygon[:len(polygon)-1]...)
516 433
434 feedback.Info("Triangulate new point cloud.")
435
436 start = time.Now()
437
517 xyz = octree.MultiPointZ(generated) 438 xyz = octree.MultiPointZ(generated)
518 439
519 start = time.Now()
520 tri, err := octree.Triangulate(xyz) 440 tri, err := octree.Triangulate(xyz)
521 if err != nil { 441 if err != nil {
522 return nil, err 442 return nil, err
523 } 443 }
524 log.Printf("Second triangulation took %v.", time.Since(start)) 444 feedback.Info("Second triangulation took %v.", time.Since(start))
525 log.Printf("Number triangles: %d.", len(tri.Triangles)/3) 445 feedback.Info("Number triangles: %d.", len(tri.Triangles)/3)
526 446 feedback.Info("Clipping triangles from new mesh.")
447
448 start = time.Now()
527 secondTin := tri.Tin() 449 secondTin := tri.Tin()
528 450 secondTin.EPSG = uint32(epsg)
529 //polygon.Buffer(1.0)
530
531 clippingPolygonBuffered.Indexify() 451 clippingPolygonBuffered.Indexify()
532 452
533 var str octree.STRTree 453 var str octree.STRTree
534 str.Build(secondTin) 454 str.Build(secondTin)
535 log.Printf("Building STR tree took %v", time.Since(start)) 455 feedback.Info("Building STR tree took %v", time.Since(start))
536 456
537 start = time.Now() 457 start = time.Now()
538 458
539 removed = str.Clip(&clippingPolygonBuffered) 459 removed = str.Clip(&clippingPolygonBuffered)
540 log.Printf("Clipping STR tree took %v.", time.Since(start)) 460 feedback.Info("Clipping STR tree took %v.", time.Since(start))
541 log.Printf("Number of triangles to clip %d.", len(removed)) 461 feedback.Info("Number of triangles to clip %d.", len(removed))
542 462
543 start = time.Now() 463 start = time.Now()
544 464
545 secondTin.EPSG = uint32(epsg) 465 feedback.Info("Building final octree index")
546 466
547 builder = octree.NewBuilder(secondTin) 467 builder = octree.NewBuilder(secondTin)
548 builder.Build(removed) 468 builder.Build(removed)
549 octreeIndex, err := builder.Bytes() 469 octreeIndex, err := builder.Bytes()
550 if err != nil { 470 if err != nil {
551 return nil, err 471 return nil, err
552 } 472 }
553 feedback.Info("Building octree took %v.", time.Since(start)) 473 feedback.Info("Building octree took %v.", time.Since(start))
474
475 start = time.Now()
554 476
555 var ( 477 var (
556 id int64 478 id int64
557 dummy uint32 479 dummy uint32
558 lat, lon float64 480 lat, lon float64
577 &hull, 499 &hull,
578 ); err != nil { 500 ); err != nil {
579 return nil, err 501 return nil, err
580 } 502 }
581 503
582 start = time.Now()
583 h := sha1.New() 504 h := sha1.New()
584 h.Write(octreeIndex) 505 h.Write(octreeIndex)
585 checksum := hex.EncodeToString(h.Sum(nil)) 506 checksum := hex.EncodeToString(h.Sum(nil))
586 _, err = tx.ExecContext(ctx, insertOctreeSQL, id, checksum, octreeIndex) 507 _, err = tx.ExecContext(ctx, insertOctreeSQL, id, checksum, octreeIndex)
587 if err != nil { 508 if err != nil {