Mercurial > gemma
comparison pkg/imports/sr.go @ 4828:39ee00d06309
Merged remove-octree-debris branch back into default.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Tue, 05 Nov 2019 14:31:22 +0100 |
parents | f4abfd0ee8ad |
children | 046a07a33b19 |
comparison
equal
deleted
inserted
replaced
4825:5eb05714353a | 4828:39ee00d06309 |
---|---|
33 "time" | 33 "time" |
34 | 34 |
35 shp "github.com/jonas-p/go-shp" | 35 shp "github.com/jonas-p/go-shp" |
36 | 36 |
37 "gemma.intevation.de/gemma/pkg/common" | 37 "gemma.intevation.de/gemma/pkg/common" |
38 "gemma.intevation.de/gemma/pkg/mesh" | |
38 "gemma.intevation.de/gemma/pkg/models" | 39 "gemma.intevation.de/gemma/pkg/models" |
39 "gemma.intevation.de/gemma/pkg/octree" | |
40 "gemma.intevation.de/gemma/pkg/wkb" | 40 "gemma.intevation.de/gemma/pkg/wkb" |
41 ) | 41 ) |
42 | 42 |
43 // SoundingResult is a Job to import sounding reults | 43 // SoundingResult is a Job to import sounding reults |
44 // from a ZIP file into the database. | 44 // from a ZIP file into the database. |
312 | 312 |
313 // LDC is cm. The data is in m. | 313 // LDC is cm. The data is in m. |
314 ldc /= 100 | 314 ldc /= 100 |
315 xform = chainTransforms( | 315 xform = chainTransforms( |
316 xform, | 316 xform, |
317 func(v octree.Vertex) octree.Vertex { | 317 func(v mesh.Vertex) mesh.Vertex { |
318 return octree.Vertex{X: v.X, Y: v.Y, Z: v.Z + ldc} | 318 return mesh.Vertex{X: v.X, Y: v.Y, Z: v.Z + ldc} |
319 }) | 319 }) |
320 m.DepthReference = depthReference | 320 m.DepthReference = depthReference |
321 } | 321 } |
322 | 322 |
323 if err := m.Validate(ctx, conn); err != nil { | 323 if err := m.Validate(ctx, conn); err != nil { |
324 return nil, common.ToError(err) | 324 return nil, common.ToError(err) |
325 } | 325 } |
326 | 326 |
327 var xyz octree.MultiPointZ | 327 var xyz mesh.MultiPointZ |
328 | 328 |
329 if z != nil { // Scanning ZIP file for *.xyz file. | 329 if z != nil { // Scanning ZIP file for *.xyz file. |
330 var xyzf *zip.File | 330 var xyzf *zip.File |
331 for _, ext := range []string{".xyz", ".txt"} { | 331 for _, ext := range []string{".xyz", ".txt"} { |
332 feedback.Info("Looking for '*%s'", ext) | 332 feedback.Info("Looking for '*%s'", ext) |
390 ctx context.Context, | 390 ctx context.Context, |
391 tx *sql.Tx, | 391 tx *sql.Tx, |
392 feedback Feedback, | 392 feedback Feedback, |
393 importID int64, | 393 importID int64, |
394 m *models.SoundingResultMeta, | 394 m *models.SoundingResultMeta, |
395 xyz octree.MultiPointZ, | 395 xyz mesh.MultiPointZ, |
396 boundary polygonSlice, | 396 boundary polygonSlice, |
397 ) (interface{}, error) { | 397 ) (interface{}, error) { |
398 | 398 |
399 if sr.singleBeam() { | 399 if sr.singleBeam() { |
400 feedback.Info("Processing as single beam scan.") | 400 feedback.Info("Processing as single beam scan.") |
427 epsg, time.Since(start)) | 427 epsg, time.Since(start)) |
428 feedback.Info("Number of reprojected points: %d", len(xyz)) | 428 feedback.Info("Number of reprojected points: %d", len(xyz)) |
429 feedback.Info("Triangulate XYZ data.") | 429 feedback.Info("Triangulate XYZ data.") |
430 | 430 |
431 start = time.Now() | 431 start = time.Now() |
432 tri, err := octree.Triangulate(xyz) | 432 tri, err := mesh.Triangulate(xyz) |
433 if err != nil { | 433 if err != nil { |
434 return nil, err | 434 return nil, err |
435 } | 435 } |
436 feedback.Info("Triangulation took %v.", time.Since(start)) | 436 feedback.Info("Triangulation took %v.", time.Since(start)) |
437 feedback.Info("Number triangles: %d.", len(tri.Triangles)/3) | 437 feedback.Info("Number triangles: %d.", len(tri.Triangles)/3) |
438 | 438 |
439 var ( | 439 var ( |
440 clippingPolygon octree.Polygon | 440 clippingPolygon mesh.Polygon |
441 clippingPolygonBuffered octree.Polygon | 441 clippingPolygonBuffered mesh.Polygon |
442 removed map[int32]struct{} | 442 removed map[int32]struct{} |
443 polygonArea float64 | 443 polygonArea float64 |
444 clippingPolygonWKB []byte | 444 clippingPolygonWKB []byte |
445 tin *octree.Tin | 445 tin *mesh.Tin |
446 ) | 446 ) |
447 | 447 |
448 if boundary == nil { | 448 if boundary == nil { |
449 feedback.Info("No boundary given. Calulate from XYZ data.") | 449 feedback.Info("No boundary given. Calulate from XYZ data.") |
450 tooLongEdge := tri.EstimateTooLong() | 450 tooLongEdge := tri.EstimateTooLong() |
451 feedback.Info("Eliminate triangles with edges longer than %.2f meters.", tooLongEdge) | 451 feedback.Info("Eliminate triangles with edges longer than %.2f meters.", tooLongEdge) |
452 | 452 |
453 var polygon octree.LineStringZ | 453 var polygon mesh.LineStringZ |
454 start = time.Now() | 454 start = time.Now() |
455 polygon, removed = tri.ConcaveHull(tooLongEdge) | 455 polygon, removed = tri.ConcaveHull(tooLongEdge) |
456 | 456 |
457 polygonArea = polygon.Area() | 457 polygonArea = polygon.Area() |
458 if polygonArea < 0.0 { // counter clockwise | 458 if polygonArea < 0.0 { // counter clockwise |
501 } | 501 } |
502 | 502 |
503 tin := tri.Tin() | 503 tin := tri.Tin() |
504 tin.EPSG = epsg | 504 tin.EPSG = epsg |
505 | 505 |
506 var str octree.STRTree | 506 var str mesh.STRTree |
507 str.Build(tin) | 507 str.Build(tin) |
508 | 508 |
509 removed = str.Clip(&clippingPolygon) | 509 removed = str.Clip(&clippingPolygon) |
510 } | 510 } |
511 | 511 |
527 feedback.Info("Build virtual DEM based on original XYZ data.") | 527 feedback.Info("Build virtual DEM based on original XYZ data.") |
528 | 528 |
529 start = time.Now() | 529 start = time.Now() |
530 | 530 |
531 tin := tri.Tin() | 531 tin := tri.Tin() |
532 var virtual octree.STRTree | 532 var virtual mesh.STRTree |
533 virtual.BuildWithout(tin, removed) | 533 virtual.BuildWithout(tin, removed) |
534 | 534 |
535 feedback.Info("Building took %v", time.Since(start)) | 535 feedback.Info("Building took %v", time.Since(start)) |
536 | 536 |
537 numPoints := int(math.Ceil(polygonArea * pointsPerSquareMeter)) | 537 numPoints := int(math.Ceil(polygonArea * pointsPerSquareMeter)) |
539 feedback.Info("Generate %d random points for an average density of ~%.2f points/m².", | 539 feedback.Info("Generate %d random points for an average density of ~%.2f points/m².", |
540 numPoints, pointsPerSquareMeter) | 540 numPoints, pointsPerSquareMeter) |
541 | 541 |
542 start = time.Now() | 542 start = time.Now() |
543 | 543 |
544 generated := make(octree.LineStringZ, 0, numPoints+clippingPolygon.NumVertices(0)) | 544 generated := make(mesh.LineStringZ, 0, numPoints+clippingPolygon.NumVertices(0)) |
545 | 545 |
546 octree.GenerateRandomVertices( | 546 mesh.GenerateRandomVertices( |
547 numPoints, | 547 numPoints, |
548 tin.Min, tin.Max, | 548 tin.Min, tin.Max, |
549 virtual.Value, | 549 virtual.Value, |
550 func(vertices []octree.Vertex) { | 550 func(vertices []mesh.Vertex) { |
551 generated = append(generated, vertices...) | 551 generated = append(generated, vertices...) |
552 }) | 552 }) |
553 | 553 |
554 feedback.Info("Generating %d points took %v.", len(generated), time.Since(start)) | 554 feedback.Info("Generating %d points took %v.", len(generated), time.Since(start)) |
555 | 555 |
560 if _, found := dupes[key]; found { | 560 if _, found := dupes[key]; found { |
561 return | 561 return |
562 } | 562 } |
563 dupes[key] = struct{}{} | 563 dupes[key] = struct{}{} |
564 if z, ok := virtual.Value(x, y); ok { | 564 if z, ok := virtual.Value(x, y); ok { |
565 generated = append(generated, octree.Vertex{X: x, Y: y, Z: z}) | 565 generated = append(generated, mesh.Vertex{X: x, Y: y, Z: z}) |
566 } | 566 } |
567 }) | 567 }) |
568 | 568 |
569 feedback.Info("Triangulate new point cloud.") | 569 feedback.Info("Triangulate new point cloud.") |
570 xyz = octree.MultiPointZ(generated) | 570 xyz = mesh.MultiPointZ(generated) |
571 start = time.Now() | 571 start = time.Now() |
572 | 572 |
573 tri, err = octree.Triangulate(xyz) | 573 tri, err = mesh.Triangulate(xyz) |
574 if err != nil { | 574 if err != nil { |
575 return nil, err | 575 return nil, err |
576 } | 576 } |
577 feedback.Info("Second triangulation took %v.", time.Since(start)) | 577 feedback.Info("Second triangulation took %v.", time.Since(start)) |
578 feedback.Info("Number triangles: %d.", len(tri.Triangles)/3) | 578 feedback.Info("Number triangles: %d.", len(tri.Triangles)/3) |
583 | 583 |
584 start = time.Now() | 584 start = time.Now() |
585 tin = tri.Tin() | 585 tin = tri.Tin() |
586 tin.EPSG = epsg | 586 tin.EPSG = epsg |
587 | 587 |
588 var str octree.STRTree | 588 var str mesh.STRTree |
589 str.Build(tin) | 589 str.Build(tin) |
590 feedback.Info("Building clipping index took %v", time.Since(start)) | 590 feedback.Info("Building clipping index took %v", time.Since(start)) |
591 | 591 |
592 start = time.Now() | 592 start = time.Now() |
593 | 593 |
595 removed = str.Clip(&clippingPolygonBuffered) | 595 removed = str.Clip(&clippingPolygonBuffered) |
596 feedback.Info("Clipping took %v.", time.Since(start)) | 596 feedback.Info("Clipping took %v.", time.Since(start)) |
597 feedback.Info("Number of triangles to clip %d.", len(removed)) | 597 feedback.Info("Number of triangles to clip %d.", len(removed)) |
598 | 598 |
599 start = time.Now() | 599 start = time.Now() |
600 final := octree.STRTree{Entries: 16} | 600 final := mesh.STRTree{Entries: 16} |
601 final.BuildWithout(tin, removed) | 601 final.BuildWithout(tin, removed) |
602 | 602 |
603 feedback.Info("Building final mesh took %v.", time.Since(start)) | 603 feedback.Info("Building final mesh took %v.", time.Since(start)) |
604 feedback.Info("Store mesh.") | 604 feedback.Info("Store mesh.") |
605 | 605 |
740 } | 740 } |
741 | 741 |
742 return &m, nil | 742 return &m, nil |
743 } | 743 } |
744 | 744 |
745 type vertexTransform func(octree.Vertex) octree.Vertex | 745 type vertexTransform func(mesh.Vertex) mesh.Vertex |
746 | 746 |
747 func identityTransform(v octree.Vertex) octree.Vertex { return v } | 747 func identityTransform(v mesh.Vertex) mesh.Vertex { return v } |
748 | 748 |
749 func negateZTransform(v octree.Vertex) octree.Vertex { | 749 func negateZTransform(v mesh.Vertex) mesh.Vertex { |
750 return octree.Vertex{X: v.X, Y: v.Y, Z: -v.Z} | 750 return mesh.Vertex{X: v.X, Y: v.Y, Z: -v.Z} |
751 } | 751 } |
752 | 752 |
753 func chainTransforms(a, b vertexTransform) vertexTransform { | 753 func chainTransforms(a, b vertexTransform) vertexTransform { |
754 return func(v octree.Vertex) octree.Vertex { return b(a(v)) } | 754 return func(v mesh.Vertex) mesh.Vertex { return b(a(v)) } |
755 } | 755 } |
756 | 756 |
757 func loadXYZReader(r io.Reader, feedback Feedback, xform vertexTransform) (octree.MultiPointZ, error) { | 757 func loadXYZReader(r io.Reader, feedback Feedback, xform vertexTransform) (mesh.MultiPointZ, error) { |
758 mpz := make(octree.MultiPointZ, 0, 250000) | 758 mpz := make(mesh.MultiPointZ, 0, 250000) |
759 s := bufio.NewScanner(r) | 759 s := bufio.NewScanner(r) |
760 | 760 |
761 warnLimiter := common.WarningLimiter{Log: feedback.Warn, MaxWarnings: 100} | 761 warnLimiter := common.WarningLimiter{Log: feedback.Warn, MaxWarnings: 100} |
762 warn := warnLimiter.Warn | 762 warn := warnLimiter.Warn |
763 defer warnLimiter.Close() | 763 defer warnLimiter.Close() |
764 | 764 |
765 for line := 1; s.Scan(); line++ { | 765 for line := 1; s.Scan(); line++ { |
766 text := s.Text() | 766 text := s.Text() |
767 var p octree.Vertex | 767 var p mesh.Vertex |
768 // fmt.Sscanf(text, "%f,%f,%f") is 4 times slower. | 768 // fmt.Sscanf(text, "%f,%f,%f") is 4 times slower. |
769 idx := strings.IndexByte(text, ',') | 769 idx := strings.IndexByte(text, ',') |
770 if idx == -1 { | 770 if idx == -1 { |
771 warn("format error in line %d", line) | 771 warn("format error in line %d", line) |
772 continue | 772 continue |
799 } | 799 } |
800 | 800 |
801 return mpz, nil | 801 return mpz, nil |
802 } | 802 } |
803 | 803 |
804 func loadXYZ(f *zip.File, feedback Feedback, xform vertexTransform) (octree.MultiPointZ, error) { | 804 func loadXYZ(f *zip.File, feedback Feedback, xform vertexTransform) (mesh.MultiPointZ, error) { |
805 r, err := f.Open() | 805 r, err := f.Open() |
806 if err != nil { | 806 if err != nil { |
807 return nil, err | 807 return nil, err |
808 } | 808 } |
809 defer r.Close() | 809 defer r.Close() |
810 return loadXYZReader(r, feedback, xform) | 810 return loadXYZReader(r, feedback, xform) |
811 } | 811 } |
812 | 812 |
813 func loadXYZFile(f string, feedback Feedback, xform vertexTransform) (octree.MultiPointZ, error) { | 813 func loadXYZFile(f string, feedback Feedback, xform vertexTransform) (mesh.MultiPointZ, error) { |
814 r, err := os.Open(f) | 814 r, err := os.Open(f) |
815 if err != nil { | 815 if err != nil { |
816 return nil, err | 816 return nil, err |
817 } | 817 } |
818 defer r.Close() | 818 defer r.Close() |
857 | 857 |
858 func generateIsos( | 858 func generateIsos( |
859 ctx context.Context, | 859 ctx context.Context, |
860 tx *sql.Tx, | 860 tx *sql.Tx, |
861 feedback Feedback, | 861 feedback Feedback, |
862 tree *octree.STRTree, | 862 tree *mesh.STRTree, |
863 id int64, | 863 id int64, |
864 ) error { | 864 ) error { |
865 | 865 |
866 heights, err := octree.LoadClassBreaks( | 866 heights, err := mesh.LoadClassBreaks( |
867 ctx, tx, | 867 ctx, tx, |
868 "morphology_classbreaks", | 868 "morphology_classbreaks", |
869 ) | 869 ) |
870 | 870 |
871 minZ, maxZ := tree.Min().Z, tree.Max().Z | 871 minZ, maxZ := tree.Min().Z, tree.Max().Z |
877 h := contourStepWidth * math.Ceil(minZ/contourStepWidth) | 877 h := contourStepWidth * math.Ceil(minZ/contourStepWidth) |
878 for ; h <= maxZ; h += contourStepWidth { | 878 for ; h <= maxZ; h += contourStepWidth { |
879 heights = append(heights, h) | 879 heights = append(heights, h) |
880 } | 880 } |
881 } else { | 881 } else { |
882 heights = octree.ExtrapolateClassBreaks(heights, minZ, maxZ) | 882 heights = mesh.ExtrapolateClassBreaks(heights, minZ, maxZ) |
883 } | 883 } |
884 | 884 |
885 /* | 885 /* |
886 for i, v := range heights { | 886 for i, v := range heights { |
887 fmt.Printf("%d %.2f\n", i, v) | 887 fmt.Printf("%d %.2f\n", i, v) |
896 | 896 |
897 func generateIsoAreas( | 897 func generateIsoAreas( |
898 ctx context.Context, | 898 ctx context.Context, |
899 tx *sql.Tx, | 899 tx *sql.Tx, |
900 feedback Feedback, | 900 feedback Feedback, |
901 tree *octree.STRTree, | 901 tree *mesh.STRTree, |
902 heights []float64, | 902 heights []float64, |
903 id int64, | 903 id int64, |
904 ) error { | 904 ) error { |
905 feedback.Info("Generate iso areas") | 905 feedback.Info("Generate iso areas") |
906 total := time.Now() | 906 total := time.Now() |
907 defer func() { | 907 defer func() { |
908 feedback.Info("Generating iso areas took %s.", | 908 feedback.Info("Generating iso areas took %s.", |
909 time.Since(total)) | 909 time.Since(total)) |
910 }() | 910 }() |
911 | 911 |
912 box := octree.Box2D{ | 912 box := mesh.Box2D{ |
913 X1: tree.Min().X, | 913 X1: tree.Min().X, |
914 Y1: tree.Min().Y, | 914 Y1: tree.Min().Y, |
915 X2: tree.Max().X, | 915 X2: tree.Max().X, |
916 Y2: tree.Max().Y, | 916 Y2: tree.Max().Y, |
917 } | 917 } |
918 | 918 |
919 raster := octree.NewRaster(box, isoCellSize) | 919 raster := mesh.NewRaster(box, isoCellSize) |
920 raster.Rasterize(tree.Value) | 920 raster.Rasterize(tree.Value) |
921 areas := raster.Trace(heights) | 921 areas := raster.Trace(heights) |
922 | 922 |
923 return storeAreas( | 923 return storeAreas( |
924 ctx, tx, feedback, | 924 ctx, tx, feedback, |