Mercurial > gemma

comparison pkg/imports/sr.go @ 3742:879c297c47e9 concave-hull

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
SR import: Use home-brew concave hull algorithm in multi-beam scan case, too.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Mon, 24 Jun 2019 12:54:25 +0200
parents ec86a7155377
children 4bb5dfa0b7e3
comparison
equal deleted inserted replaced
3733:ec86a7155377 3742:879c297c47e9
304 defer tx.Rollback() 304 defer tx.Rollback()
305 305
306 var summary interface{} 306 var summary interface{}
307 307
308 if sr.isSingleBeam() { 308 if sr.isSingleBeam() {
309 summary, err = sr.singleBeamScan( 309 summary, err = sr.processScan(
310 ctx, 310 ctx,
311 tx, 311 tx,
312 feedback, 312 feedback,
313 true,
313 importID, 314 importID,
314 m, 315 m,
315 xyz, 316 xyz,
316 boundary, 317 boundary,
317 ) 318 )
318 } else { 319 } else {
319 summary, err = sr.multiBeamScan( 320 summary, err = sr.processScan(
320 ctx, 321 ctx,
321 tx, 322 tx,
322 feedback, 323 feedback,
324 false,
323 importID, 325 importID,
324 m, 326 m,
325 xyz, 327 xyz,
326 boundary, 328 boundary,
327 ) 329 )
340 "Storing sounding result was successful after %v.", time.Since(start)) 342 "Storing sounding result was successful after %v.", time.Since(start))
341 343
342 return summary, nil 344 return summary, nil
343 } 345 }
344 346
345 func (sr *SoundingResult) singleBeamScan( 347 func (sr *SoundingResult) processScan(
346 ctx context.Context, 348 ctx context.Context,
347 tx *sql.Tx, 349 tx *sql.Tx,
348 feedback Feedback, 350 feedback Feedback,
351 singleBeam bool,
349 importID int64, 352 importID int64,
350 m *models.SoundingResultMeta, 353 m *models.SoundingResultMeta,
351 xyz octree.MultiPointZ, 354 xyz octree.MultiPointZ,
352 boundary polygonSlice, 355 boundary polygonSlice,
353 ) (interface{}, error) { 356 ) (interface{}, error) {
354 357
355 feedback.Info("Processing as single beam scan.") 358 if singleBeam {
359 feedback.Info("Processing as single beam scan.")
360 } else {
361 feedback.Info("Processing as multi beam scan.")
362 }
363
356 feedback.Info("Reproject XYZ data.") 364 feedback.Info("Reproject XYZ data.")
357 365
358 start := time.Now() 366 start := time.Now()
359 367
360 xyzWKB := xyz.AsWKB() 368 xyzWKB := xyz.AsWKB()
361 var reproj []byte 369 var reproj []byte
362 var epsg uint 370 var epsg uint32
363 371
364 if err := tx.QueryRowContext( 372 if err := tx.QueryRowContext(
365 ctx, 373 ctx,
366 reprojectPointsSingleBeamSQL, 374 reprojectPointsSingleBeamSQL,
367 xyzWKB, 375 xyzWKB,
391 clippingPolygon octree.Polygon 399 clippingPolygon octree.Polygon
392 clippingPolygonBuffered octree.Polygon 400 clippingPolygonBuffered octree.Polygon
393 removed map[int32]struct{} 401 removed map[int32]struct{}
394 polygonArea float64 402 polygonArea float64
395 clippingPolygonWKB []byte 403 clippingPolygonWKB []byte
404 tin *octree.Tin
396 ) 405 )
397 406
398 if boundary == nil { 407 if boundary == nil {
399 feedback.Info("No boundary given. Calulate from XYZ data.") 408 feedback.Info("No boundary given. Calulate from XYZ data.")
400 tooLongEdge := tri.EstimateTooLong() 409 tooLongEdge := tri.EstimateTooLong()
449 if err := clippingPolygonBuffered.FromWKB(buffered); err != nil { 458 if err := clippingPolygonBuffered.FromWKB(buffered); err != nil {
450 return nil, err 459 return nil, err
451 } 460 }
452 461
453 tin := tri.Tin() 462 tin := tri.Tin()
454 tin.EPSG = uint32(epsg) 463 tin.EPSG = epsg
455 464
456 var str octree.STRTree 465 var str octree.STRTree
457 str.Build(tin) 466 str.Build(tin)
458 467
459 removed = str.Clip(&clippingPolygon) 468 removed = str.Clip(&clippingPolygon)
460 } 469 }
461 470
462 // Build the first mesh to generate random points on. 471 if singleBeam {
463 472
464 feedback.Info("Build virtual DEM based on original XYZ data.") 473 // Build the first mesh to generate random points on.
474
475 feedback.Info("Build virtual DEM based on original XYZ data.")
476
477 start = time.Now()
478
479 var tree *octree.Tree
480 {
481 builder := octree.NewBuilder(tri.Tin())
482 builder.Build(removed)
483 tree = builder.Tree()
484 }
485
486 feedback.Info("Building took %v", time.Since(start))
487
488 feedback.Info("Boundary area: %.2fm²", polygonArea)
489
490 numPoints := int(math.Ceil(polygonArea * pointsPerSquareMeter))
491
492 feedback.Info("Generate %d random points for an average density of ~%d points/m².",
493 numPoints, pointsPerSquareMeter)
494
495 start = time.Now()
496
497 generated := make(octree.LineStringZ, 0, numPoints+clippingPolygon.NumVertices(0))
498
499 tree.GenerateRandomVertices(numPoints, func(vertices []octree.Vertex) {
500 generated = append(generated, vertices...)
501 })
502
503 feedback.Info("Generating %d points took %v.", len(generated), time.Since(start))
504
505 // Add the boundary to new point cloud.
506 dupes := map[[2]float64]struct{}{}
507 clippingPolygon.Vertices(0, func(x, y float64) {
508 key := [2]float64{x, y}
509 if _, found := dupes[key]; found {
510 return
511 }
512 dupes[key] = struct{}{}
513 if z, ok := tree.Value(x, y); ok {
514 generated = append(generated, octree.Vertex{X: x, Y: y, Z: z})
515 }
516 })
517
518 feedback.Info("Triangulate new point cloud.")
519 xyz = octree.MultiPointZ(generated)
520 start = time.Now()
521
522 tri, err = octree.Triangulate(xyz)
523 if err != nil {
524 return nil, err
525 }
526 feedback.Info("Second triangulation took %v.", time.Since(start))
527 feedback.Info("Number triangles: %d.", len(tri.Triangles)/3)
528 feedback.Info("Clipping triangles from new mesh.")
529
530 } else { // multi beam
531 // Nothing special
532 }
465 533
466 start = time.Now() 534 start = time.Now()
467 535 tin = tri.Tin()
468 var tree *octree.Tree 536 tin.EPSG = epsg
469 {
470 builder := octree.NewBuilder(tri.Tin())
471 builder.Build(removed)
472 tree = builder.Tree()
473 }
474
475 feedback.Info("Building took %v", time.Since(start))
476
477 feedback.Info("Boundary area: %.2fm²", polygonArea)
478
479 numPoints := int(math.Ceil(polygonArea * pointsPerSquareMeter))
480
481 feedback.Info("Generate %d random points for an average density of ~%d points/m².",
482 numPoints, pointsPerSquareMeter)
483
484 start = time.Now()
485
486 generated := make(octree.LineStringZ, 0, numPoints+clippingPolygon.NumVertices(0))
487
488 tree.GenerateRandomVertices(numPoints, func(vertices []octree.Vertex) {
489 generated = append(generated, vertices...)
490 })
491
492 feedback.Info("Generating %d points took %v.", len(generated), time.Since(start))
493
494 // Add the boundary to new point cloud.
495 dupes := map[[2]float64]struct{}{}
496 clippingPolygon.Vertices(0, func(x, y float64) {
497 key := [2]float64{x, y}
498 if _, found := dupes[key]; found {
499 return
500 }
501 dupes[key] = struct{}{}
502 if z, ok := tree.Value(x, y); ok {
503 generated = append(generated, octree.Vertex{X: x, Y: y, Z: z})
504 }
505 })
506
507 feedback.Info("Triangulate new point cloud.")
508
509 start = time.Now()
510
511 xyz = octree.MultiPointZ(generated)
512
513 tri, err = octree.Triangulate(xyz)
514 if err != nil {
515 return nil, err
516 }
517 feedback.Info("Second triangulation took %v.", time.Since(start))
518 feedback.Info("Number triangles: %d.", len(tri.Triangles)/3)
519 feedback.Info("Clipping triangles from new mesh.")
520
521 start = time.Now()
522 tin := tri.Tin()
523 tin.EPSG = uint32(epsg)
524 537
525 var str octree.STRTree 538 var str octree.STRTree
526 str.Build(tin) 539 str.Build(tin)
527 feedback.Info("Building STR tree took %v", time.Since(start)) 540 feedback.Info("Building STR tree took %v", time.Since(start))
528 541
530 543
531 clippingPolygonBuffered.Indexify() 544 clippingPolygonBuffered.Indexify()
532 removed = str.Clip(&clippingPolygonBuffered) 545 removed = str.Clip(&clippingPolygonBuffered)
533 feedback.Info("Clipping STR tree took %v.", time.Since(start)) 546 feedback.Info("Clipping STR tree took %v.", time.Since(start))
534 feedback.Info("Number of triangles to clip %d.", len(removed)) 547 feedback.Info("Number of triangles to clip %d.", len(removed))
535
536 start = time.Now()
537 548
538 feedback.Info("Build final octree index") 549 feedback.Info("Build final octree index")
539 550
540 builder := octree.NewBuilder(tin) 551 builder := octree.NewBuilder(tin)
541 builder.Build(removed) 552 builder.Build(removed)
617 } 628 }
618 629
619 return &summary, nil 630 return &summary, nil
620 } 631 }
621 632
622 func (sr *SoundingResult) multiBeamScan(
623 ctx context.Context,
624 tx *sql.Tx,
625 feedback Feedback,
626 importID int64,
627 m *models.SoundingResultMeta,
628 xyz octree.MultiPointZ,
629 boundary polygonSlice,
630 ) (interface{}, error) {
631
632 feedback.Info("Processing as multi beam scan.")
633 var (
634 id int64
635 epsg uint32
636 lat, lon float64
637 )
638 start := time.Now()
639
640 var hull []byte
641
642 xyzWKB := xyz.AsWKB()
643
644 err := tx.QueryRowContext(
645 ctx,
646 insertHullSQL,
647 m.Bottleneck,
648 m.Date.Time,
649 m.DepthReference,
650 xyzWKB,
651 boundary.asWKB(),
652 m.EPSG,
653 ).Scan(
654 &id,
655 &lat,
656 &lon,
657 &epsg,
658 &hull,
659 )
660 xyz, boundary = nil, nil // not need from now on.
661 feedback.Info("Calculating hull took %s.", time.Since(start))
662 switch {
663 case err == sql.ErrNoRows:
664 return nil, fmt.Errorf(
665 "No matching bottleneck of given name or time available: %v", err)
666 case err != nil:
667 return nil, err
668 }
669 feedback.Info("Best suited UTM EPSG: %d", epsg)
670
671 start = time.Now()
672
673 var clippingPolygon octree.Polygon
674
675 if err := clippingPolygon.FromWKB(hull); err != nil {
676 return nil, err
677 }
678 clippingPolygon.Indexify()
679 feedback.Info("Building clipping polygon took %v.", time.Since(start))
680
681 start = time.Now()
682
683 var reproj []byte
684
685 if err = tx.QueryRowContext(
686 ctx,
687 reprojectPointsSQL,
688 xyzWKB,
689 m.EPSG,
690 epsg,
691 ).Scan(&reproj); err != nil {
692 return nil, err
693 }
694
695 if err := xyz.FromWKB(reproj); err != nil {
696 return nil, err
697 }
698
699 feedback.Info("Reprojecting points took %v.", time.Since(start))
700 feedback.Info("Number of reprojected points: %d", len(xyz))
701
702 start = time.Now()
703
704 tri, err := octree.Triangulate(xyz)
705 if err != nil {
706 return nil, err
707 }
708 feedback.Info("Triangulation took %v.", time.Since(start))
709
710 start = time.Now()
711
712 tin := tri.Tin()
713
714 var str octree.STRTree
715 str.Build(tin)
716 feedback.Info("Building STR tree took %v", time.Since(start))
717
718 start = time.Now()
719
720 removed := str.Clip(&clippingPolygon)
721 feedback.Info("Clipping STR tree took %v.", time.Since(start))
722 feedback.Info("Number of triangles to clip %d.", len(removed))
723
724 start = time.Now()
725
726 tin.EPSG = epsg
727
728 builder := octree.NewBuilder(tin)
729 builder.Build(removed)
730 octreeIndex, err := builder.Bytes()
731 if err != nil {
732 return nil, err
733 }
734 feedback.Info("Building octree took %v.", time.Since(start))
735
736 start = time.Now()
737 h := sha1.New()
738 h.Write(octreeIndex)
739 checksum := hex.EncodeToString(h.Sum(nil))
740 _, err = tx.ExecContext(ctx, insertOctreeSQL, id, checksum, octreeIndex)
741 if err != nil {
742 return nil, err
743 }
744 feedback.Info("Storing octree index took %s.", time.Since(start))
745
746 tree := builder.Tree()
747
748 start = time.Now()
749 err = generateContours(ctx, tx, tree, id)
750 if err != nil {
751 return nil, err
752 }
753 feedback.Info("Generating and storing contour lines took %s.",
754 time.Since(start))
755
756 // Store for potential later removal.
757 if err = track(ctx, tx, importID, "waterway.sounding_results", id); err != nil {
758 return nil, err
759 }
760
761 summary := struct {
762 Bottleneck string `json:"bottleneck"`
763 Date models.Date `json:"date"`
764 Lat float64 `json:"lat"`
765 Lon float64 `json:"lon"`
766 }{
767 Bottleneck: m.Bottleneck,
768 Date: m.Date,
769 Lat: lat,
770 Lon: lon,
771 }
772
773 return &summary, nil
774 }
775
776 // CleanUp removes the folder containing the ZIP file with the 633 // CleanUp removes the folder containing the ZIP file with the
777 // the sounding result import. 634 // the sounding result import.
778 func (sr *SoundingResult) CleanUp() error { 635 func (sr *SoundingResult) CleanUp() error {
779 return os.RemoveAll(sr.Dir) 636 return os.RemoveAll(sr.Dir)
780 } 637 }