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 |