Mercurial > gemma
comparison pkg/imports/sr.go @ 4751:dbf07d0c364e
SR import: If the point density is greater than 0.2 values per meter² the interpolating step is omitted.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Fri, 18 Oct 2019 17:07:24 +0200 |
parents | 56bd9ba0354c |
children | 64979fec89a7 |
comparison
equal
deleted
inserted
replaced
4750:67a5de15490b | 4751:dbf07d0c364e |
---|---|
66 } | 66 } |
67 | 67 |
68 const ( | 68 const ( |
69 contourStepWidth = 0.1 | 69 contourStepWidth = 0.1 |
70 contourTolerance = 0.1 | 70 contourTolerance = 0.1 |
71 ) | |
72 | |
73 const ( | |
74 // multiBeamThreshold is the number of points m² which | |
75 // is assumed that greater values are indicators for | |
76 // an already interpolated point cloud. | |
77 multiBeamThreshold = 1.0 / 5.0 | |
71 ) | 78 ) |
72 | 79 |
73 const ( | 80 const ( |
74 // pointsPerSquareMeter is the average number of points | 81 // pointsPerSquareMeter is the average number of points |
75 // when generating a artifical height model for single beam scans. | 82 // when generating a artifical height model for single beam scans. |
486 removed = str.Clip(&clippingPolygon) | 493 removed = str.Clip(&clippingPolygon) |
487 } | 494 } |
488 | 495 |
489 if sr.singleBeam() { | 496 if sr.singleBeam() { |
490 | 497 |
498 origDensity := float64(len(xyz)) / polygonArea | |
499 | |
500 feedback.Info("Boundary area: %.2fm²", polygonArea) | |
501 feedback.Info("Original point density: %.2f points/m²", origDensity) | |
502 | |
503 if origDensity > multiBeamThreshold { | |
504 feedback.Warn("The density is greater than %.2f points/m².") | |
505 feedback.Warn("It is assumed that the data is already interpolated.") | |
506 goto multibeam | |
507 } | |
508 | |
491 // Build the first mesh to generate random points on. | 509 // Build the first mesh to generate random points on. |
492 | 510 |
493 feedback.Info("Build virtual DEM based on original XYZ data.") | 511 feedback.Info("Build virtual DEM based on original XYZ data.") |
494 | 512 |
495 start = time.Now() | 513 start = time.Now() |
497 tin := tri.Tin() | 515 tin := tri.Tin() |
498 var virtual octree.STRTree | 516 var virtual octree.STRTree |
499 virtual.BuildWithout(tin, removed) | 517 virtual.BuildWithout(tin, removed) |
500 | 518 |
501 feedback.Info("Building took %v", time.Since(start)) | 519 feedback.Info("Building took %v", time.Since(start)) |
502 | |
503 feedback.Info("Boundary area: %.2fm²", polygonArea) | |
504 | 520 |
505 numPoints := int(math.Ceil(polygonArea * pointsPerSquareMeter)) | 521 numPoints := int(math.Ceil(polygonArea * pointsPerSquareMeter)) |
506 | 522 |
507 feedback.Info("Generate %d random points for an average density of ~%d points/m².", | 523 feedback.Info("Generate %d random points for an average density of ~%d points/m².", |
508 numPoints, pointsPerSquareMeter) | 524 numPoints, pointsPerSquareMeter) |
544 } | 560 } |
545 feedback.Info("Second triangulation took %v.", time.Since(start)) | 561 feedback.Info("Second triangulation took %v.", time.Since(start)) |
546 feedback.Info("Number triangles: %d.", len(tri.Triangles)/3) | 562 feedback.Info("Number triangles: %d.", len(tri.Triangles)/3) |
547 feedback.Info("Clipping triangles from new mesh.") | 563 feedback.Info("Clipping triangles from new mesh.") |
548 } | 564 } |
565 | |
566 multibeam: | |
549 | 567 |
550 start = time.Now() | 568 start = time.Now() |
551 tin = tri.Tin() | 569 tin = tri.Tin() |
552 tin.EPSG = epsg | 570 tin.EPSG = epsg |
553 | 571 |