Mercurial > gemma

comparison pkg/imports/sr.go @ 3658:1c3df921361d single-beam

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Handle th case that a boundary polygon is uploaded along side with the single beam scan.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Fri, 14 Jun 2019 11:59:14 +0200
parents a6c671abbc35
children 66f2cb789905
comparison
equal deleted inserted replaced
3655:a6c671abbc35 3658:1c3df921361d
134 best_utm(area), 134 best_utm(area),
135 ST_AsBinary(ST_Transform(area::geometry, best_utm(area))) 135 ST_AsBinary(ST_Transform(area::geometry, best_utm(area)))
136 ` 136 `
137 137
138 reprojectPointsSQL = ` 138 reprojectPointsSQL = `
139 SELECT ST_AsBinary(ST_Transform(ST_GeomFromWKB($1, $2::integer), $3::integer)) 139 SELECT ST_AsBinary(ST_Transform(ST_GeomFromWKB($1, $2::integer), $3::integer))`
140 ` 140
141 reprojectPointsBufferedSQL = `
142 SELECT
143 ST_AsBinary(ST_Transform(ST_GeomFromWKB($1, $2::integer), $3::integer)),
144 ST_AsBinary(ST_Buffer(ST_Transform(ST_GeomFromWKB($1, $2::integer), $3::integer), 0.1)),
145 ST_Area(ST_Transform(ST_GeomFromWKB($1, $2::integer), $3::integer))`
146
141 insertOctreeSQL = ` 147 insertOctreeSQL = `
142 UPDATE waterway.sounding_results SET 148 UPDATE waterway.sounding_results SET
143 octree_checksum = $2, octree_index = $3 149 octree_checksum = $2, octree_index = $3
144 WHERE id = $1` 150 WHERE id = $1`
145 151
336 xyz octree.MultiPointZ, 342 xyz octree.MultiPointZ,
337 boundary polygonSlice, 343 boundary polygonSlice,
338 ) (interface{}, error) { 344 ) (interface{}, error) {
339 345
340 feedback.Info("Processing as single beam scan.") 346 feedback.Info("Processing as single beam scan.")
347 feedback.Info("Reproject XYZ data.")
341 348
342 start := time.Now() 349 start := time.Now()
343 350
344 xyzWKB := xyz.AsWKB() 351 xyzWKB := xyz.AsWKB()
345 var reproj []byte 352 var reproj []byte
359 } 366 }
360 367
361 feedback.Info("Reprojecting points to EPSG %d took %v.", 368 feedback.Info("Reprojecting points to EPSG %d took %v.",
362 epsg, time.Since(start)) 369 epsg, time.Since(start))
363 feedback.Info("Number of reprojected points: %d", len(xyz)) 370 feedback.Info("Number of reprojected points: %d", len(xyz))
371 feedback.Info("Triangulate XYZ data.")
364 372
365 start = time.Now() 373 start = time.Now()
366 tri, err := octree.Triangulate(xyz) 374 tri, err := octree.Triangulate(xyz)
367 if err != nil { 375 if err != nil {
368 return nil, err 376 return nil, err
369 } 377 }
370 feedback.Info("Triangulation took %v.", time.Since(start)) 378 feedback.Info("Triangulation took %v.", time.Since(start))
371 feedback.Info("Number triangles: %d.", len(tri.Triangles)/3) 379 feedback.Info("Number triangles: %d.", len(tri.Triangles)/3)
372 380
373 var clippingPolygon octree.Polygon 381 var (
382 clippingPolygon octree.Polygon
383 clippingPolygonBuffered octree.Polygon
384 removed map[int32]struct{}
385 polygonArea float64
386 clippingPolygonWKB []byte
387 )
374 388
375 if boundary == nil { 389 if boundary == nil {
376 feedback.Info("No boundary given. Calulating from XYZ data.") 390 feedback.Info("No boundary given. Calulate from XYZ data.")
377 feedback.Info("Eliminate triangles with edges longer than %.2f meters.", tooLongEdge) 391 feedback.Info("Eliminate triangles with edges longer than %.2f meters.", tooLongEdge)
378 392
379 polygon, removed := tri.ConcaveHull(tooLongEdge) 393 var polygon octree.LineStringZ
380 394 start = time.Now()
381 polygonArea := polygon.Area() 395 polygon, removed = tri.ConcaveHull(tooLongEdge)
396
397 polygonArea = polygon.Area()
382 if polygonArea < 0.0 { // counter clockwise 398 if polygonArea < 0.0 { // counter clockwise
383 polygonArea = -polygonArea 399 polygonArea = -polygonArea
384 polygon.Reverse() 400 polygon.Reverse()
385 } 401 }
386 402
387 clippingPolygon.FromLineStringZ(polygon) 403 clippingPolygon.FromLineStringZ(polygon)
388 404
389 var repaired, buffered []byte 405 var buffered []byte
390
391 if err := tx.QueryRowContext( 406 if err := tx.QueryRowContext(
392 ctx, 407 ctx,
393 repairBoundarySQL, 408 repairBoundarySQL,
394 clippingPolygon.AsWKB(), 409 clippingPolygon.AsWKB(),
395 epsg, 410 epsg,
396 ).Scan(&repaired, &buffered); err != nil { 411 ).Scan(&clippingPolygonWKB, &buffered); err != nil {
397 return nil, err 412 return nil, err
398 } 413 }
399 414
400 if err := clippingPolygon.FromWKB(repaired); err != nil { 415 if err := clippingPolygon.FromWKB(clippingPolygonWKB); err != nil {
401 return nil, err 416 return nil, err
402 } 417 }
403 var clippingPolygonBuffered octree.Polygon
404 if err := clippingPolygonBuffered.FromWKB(buffered); err != nil { 418 if err := clippingPolygonBuffered.FromWKB(buffered); err != nil {
405 return nil, err 419 return nil, err
406 } 420 }
407 421
408 firstTin := tri.Tin() 422 feedback.Info("Calculating took %v.", time.Since(start))
409 builder := octree.NewBuilder(firstTin) 423
410 builder.Build(removed) 424 } else { // has Boundary
411 425 feedback.Info("Using uploaded boundary polygon.")
412 firstTree := builder.Tree() 426 var buffered []byte
413
414 feedback.Info("Boundary area: %.2fm²", polygonArea)
415
416 numPoints := int(math.Ceil(polygonArea * pointsPerSquareMeter))
417
418 feedback.Info("Generating %d random points for an average density of ~%d points/m².",
419 numPoints, pointsPerSquareMeter)
420
421 start := time.Now()
422
423 generated := make(octree.LineStringZ, 0, numPoints+len(polygon)-1)
424
425 firstTree.GenerateRandomVertices(numPoints, func(vertices []octree.Vertex) {
426 generated = append(generated, vertices...)
427 })
428
429 feedback.Info("Generating %d points took %v.", len(generated), time.Since(start))
430
431 // Add the boundary to new point cloud.
432 generated = append(generated, polygon[:len(polygon)-1]...)
433
434 feedback.Info("Triangulate new point cloud.")
435
436 start = time.Now()
437
438 xyz = octree.MultiPointZ(generated)
439
440 tri, err := octree.Triangulate(xyz)
441 if err != nil {
442 return nil, err
443 }
444 feedback.Info("Second triangulation took %v.", time.Since(start))
445 feedback.Info("Number triangles: %d.", len(tri.Triangles)/3)
446 feedback.Info("Clipping triangles from new mesh.")
447
448 start = time.Now()
449 secondTin := tri.Tin()
450 secondTin.EPSG = uint32(epsg)
451 clippingPolygonBuffered.Indexify()
452
453 var str octree.STRTree
454 str.Build(secondTin)
455 feedback.Info("Building STR tree took %v", time.Since(start))
456
457 start = time.Now()
458
459 removed = str.Clip(&clippingPolygonBuffered)
460 feedback.Info("Clipping STR tree took %v.", time.Since(start))
461 feedback.Info("Number of triangles to clip %d.", len(removed))
462
463 start = time.Now()
464
465 feedback.Info("Building final octree index")
466
467 builder = octree.NewBuilder(secondTin)
468 builder.Build(removed)
469 octreeIndex, err := builder.Bytes()
470 if err != nil {
471 return nil, err
472 }
473 feedback.Info("Building octree took %v.", time.Since(start))
474
475 start = time.Now()
476
477 var (
478 id int64
479 dummy uint32
480 lat, lon float64
481 )
482
483 var hull []byte
484
485 if err := tx.QueryRowContext(
486 ctx,
487 insertHullSQL,
488 m.Bottleneck,
489 m.Date.Time,
490 m.DepthReference,
491 nil,
492 repaired,
493 epsg,
494 ).Scan(
495 &id,
496 &lat,
497 &lon,
498 &dummy,
499 &hull,
500 ); err != nil {
501 return nil, err
502 }
503
504 h := sha1.New()
505 h.Write(octreeIndex)
506 checksum := hex.EncodeToString(h.Sum(nil))
507 _, err = tx.ExecContext(ctx, insertOctreeSQL, id, checksum, octreeIndex)
508 if err != nil {
509 return nil, err
510 }
511 feedback.Info("Storing octree index took %s.", time.Since(start))
512
513 tree := builder.Tree()
514
515 start = time.Now()
516 err = generateContours(ctx, tx, tree, id)
517 if err != nil {
518 return nil, err
519 }
520 feedback.Info("Generating and storing contour lines took %s.",
521 time.Since(start))
522
523 // Store for potential later removal.
524 if err = track(ctx, tx, importID, "waterway.sounding_results", id); err != nil {
525 return nil, err
526 }
527
528 summary := struct {
529 Bottleneck string `json:"bottleneck"`
530 Date models.Date `json:"date"`
531 Lat float64 `json:"lat"`
532 Lon float64 `json:"lon"`
533 }{
534 Bottleneck: m.Bottleneck,
535 Date: m.Date,
536 Lat: lat,
537 Lon: lon,
538 }
539
540 return &summary, nil
541
542 } else {
543 var hull []byte
544 if err = tx.QueryRowContext( 427 if err = tx.QueryRowContext(
545 ctx, 428 ctx,
546 reprojectPointsSQL, 429 reprojectPointsBufferedSQL,
547 boundary.asWKB(), 430 boundary.asWKB(),
548 m.EPSG, 431 m.EPSG,
549 epsg, 432 epsg,
550 ).Scan(&hull); err != nil { 433 ).Scan(&clippingPolygonWKB, &buffered, &polygonArea); err != nil {
551 return nil, err 434 return nil, err
552 } 435 }
553 if err := clippingPolygon.FromWKB(hull); err != nil { 436 if err := clippingPolygon.FromWKB(clippingPolygonWKB); err != nil {
554 return nil, err 437 return nil, err
555 } 438 }
556 clippingPolygon.Indexify() 439 if err := clippingPolygonBuffered.FromWKB(buffered); err != nil {
557 } 440 return nil, err
558 441 }
559 // TODO: Implement me! 442
560 return nil, errors.New("Not implemented, yet!") 443 tin := tri.Tin()
444 tin.EPSG = uint32(epsg)
445
446 var str octree.STRTree
447 str.Build(tin)
448
449 removed = str.Clip(&clippingPolygon)
450 }
451
452 // Build the first mesh to generate random points on.
453
454 feedback.Info("Build virtual DEM based on original XYZ data.")
455
456 start = time.Now()
457
458 var tree *octree.Tree
459 {
460 builder := octree.NewBuilder(tri.Tin())
461 builder.Build(removed)
462 tree = builder.Tree()
463 }
464
465 feedback.Info("Building took %v", time.Since(start))
466
467 feedback.Info("Boundary area: %.2fm²", polygonArea)
468
469 numPoints := int(math.Ceil(polygonArea * pointsPerSquareMeter))
470
471 feedback.Info("Generate %d random points for an average density of ~%d points/m².",
472 numPoints, pointsPerSquareMeter)
473
474 start = time.Now()
475
476 generated := make(octree.LineStringZ, 0, numPoints+clippingPolygon.NumVertices(0))
477
478 tree.GenerateRandomVertices(numPoints, func(vertices []octree.Vertex) {
479 generated = append(generated, vertices...)
480 })
481
482 feedback.Info("Generating %d points took %v.", len(generated), time.Since(start))
483
484 // Add the boundary to new point cloud.
485 dupes := map[[2]float64]struct{}{}
486 clippingPolygon.Vertices(0, func(x, y float64) {
487 key := [2]float64{x, y}
488 if _, found := dupes[key]; found {
489 return
490 }
491 dupes[key] = struct{}{}
492 if z, ok := tree.Value(x, y); ok {
493 generated = append(generated, octree.Vertex{X: x, Y: y, Z: z})
494 }
495 })
496
497 feedback.Info("Triangulate new point cloud.")
498
499 start = time.Now()
500
501 xyz = octree.MultiPointZ(generated)
502
503 tri, err = octree.Triangulate(xyz)
504 if err != nil {
505 return nil, err
506 }
507 feedback.Info("Second triangulation took %v.", time.Since(start))
508 feedback.Info("Number triangles: %d.", len(tri.Triangles)/3)
509 feedback.Info("Clipping triangles from new mesh.")
510
511 start = time.Now()
512 tin := tri.Tin()
513 tin.EPSG = uint32(epsg)
514
515 var str octree.STRTree
516 str.Build(tin)
517 feedback.Info("Building STR tree took %v", time.Since(start))
518
519 start = time.Now()
520
521 clippingPolygonBuffered.Indexify()
522 removed = str.Clip(&clippingPolygonBuffered)
523 feedback.Info("Clipping STR tree took %v.", time.Since(start))
524 feedback.Info("Number of triangles to clip %d.", len(removed))
525
526 start = time.Now()
527
528 feedback.Info("Build final octree index")
529
530 builder := octree.NewBuilder(tin)
531 builder.Build(removed)
532 octreeIndex, err := builder.Bytes()
533 if err != nil {
534 return nil, err
535 }
536 feedback.Info("Building octree took %v.", time.Since(start))
537 feedback.Info("Store octree.")
538
539 start = time.Now()
540
541 var (
542 id int64
543 dummy uint32
544 lat, lon float64
545 )
546
547 var hull []byte
548
549 if err := tx.QueryRowContext(
550 ctx,
551 insertHullSQL,
552 m.Bottleneck,
553 m.Date.Time,
554 m.DepthReference,
555 nil,
556 clippingPolygonWKB,
557 epsg,
558 ).Scan(
559 &id,
560 &lat,
561 &lon,
562 &dummy,
563 &hull,
564 ); err != nil {
565 return nil, err
566 }
567
568 h := sha1.New()
569 h.Write(octreeIndex)
570 checksum := hex.EncodeToString(h.Sum(nil))
571 _, err = tx.ExecContext(ctx, insertOctreeSQL, id, checksum, octreeIndex)
572 if err != nil {
573 return nil, err
574 }
575 feedback.Info("Storing octree index took %s.", time.Since(start))
576 feedback.Info("Generate contour lines")
577
578 start = time.Now()
579 err = generateContours(ctx, tx, builder.Tree(), id)
580 if err != nil {
581 return nil, err
582 }
583 feedback.Info("Generating contour lines took %s.",
584 time.Since(start))
585
586 // Store for potential later removal.
587 if err = track(ctx, tx, importID, "waterway.sounding_results", id); err != nil {
588 return nil, err
589 }
590
591 summary := struct {
592 Bottleneck string `json:"bottleneck"`
593 Date models.Date `json:"date"`
594 Lat float64 `json:"lat"`
595 Lon float64 `json:"lon"`
596 }{
597 Bottleneck: m.Bottleneck,
598 Date: m.Date,
599 Lat: lat,
600 Lon: lon,
601 }
602
603 return &summary, nil
561 } 604 }
562 605
563 func (sr *SoundingResult) multiBeamScan( 606 func (sr *SoundingResult) multiBeamScan(
564 ctx context.Context, 607 ctx context.Context,
565 tx *sql.Tx, 608 tx *sql.Tx,