Mercurial > gemma

comparison pkg/imports/sr.go @ 5416:31b0e865e7a0 marking-single-beam

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Stored marking points in database.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Wed, 07 Jul 2021 14:57:30 +0200
parents 4ad68ab239b7
children 2d294ad81241
comparison
equal deleted inserted replaced
5415:4ad68ab239b7 5416:31b0e865e7a0
22 "database/sql" 22 "database/sql"
23 "encoding/hex" 23 "encoding/hex"
24 "errors" 24 "errors"
25 "fmt" 25 "fmt"
26 "io" 26 "io"
27 "log"
27 "math" 28 "math"
28 "os" 29 "os"
29 "path" 30 "path"
30 "path/filepath" 31 "path/filepath"
31 "strconv" 32 "strconv"
196 4326 197 4326
197 ) 198 )
198 FROM waterway.sounding_results sr 199 FROM waterway.sounding_results sr
199 WHERE id = $1 200 WHERE id = $1
200 ` 201 `
202 insertMarkingPointsSQL = `
203 INSERT INTO waterway.sounding_results_marking_points (
204 sounding_result_id,
205 height,
206 points
207 )
208 SELECT
209 $1,
210 $2,
211 ST_Transform(ST_Force2D(ST_GeomFromWKB($4, $3::integer)), 4326)
212 FROM waterway.sounding_results sr
213 WHERE id = $1
214 `
201 215
202 selectGaugeLDCSQL = ` 216 selectGaugeLDCSQL = `
203 SELECT 217 SELECT
204 grwl.value, 218 grwl.value,
205 grwl.depth_reference 219 grwl.depth_reference
429 443
430 feedback.Info("Reproject XYZ data.") 444 feedback.Info("Reproject XYZ data.")
431 445
432 start := time.Now() 446 start := time.Now()
433 447
434 xyzWKB := xyz.AsWKB()
435 var reproj []byte 448 var reproj []byte
436 var epsg uint32 449 var epsg uint32
437 450
438 if err := tx.QueryRowContext( 451 if err := tx.QueryRowContext(
439 ctx, 452 ctx,
440 reprojectPointsSingleBeamSQL, 453 reprojectPointsSingleBeamSQL,
441 xyzWKB, 454 xyz.AsWKB(),
442 m.EPSG, 455 m.EPSG,
443 ).Scan(&reproj, &epsg); err != nil { 456 ).Scan(&reproj, &epsg); err != nil {
444 return nil, err 457 return nil, err
445 } 458 }
446 459
681 _, err = tx.ExecContext(ctx, insertMeshSQL, id, checksum, index) 694 _, err = tx.ExecContext(ctx, insertMeshSQL, id, checksum, index)
682 if err != nil { 695 if err != nil {
683 return nil, err 696 return nil, err
684 } 697 }
685 feedback.Info("Storing mesh index took %s.", time.Since(start)) 698 feedback.Info("Storing mesh index took %s.", time.Since(start))
686 err = generateIsos(ctx, tx, feedback, &final, id) 699 if err := generateIsos(ctx, tx, feedback, &final, id); err != nil {
687 if err != nil {
688 return nil, err 700 return nil, err
689 } 701 }
690 } else { // SurveyTypeMarking 702 } else { // SurveyTypeMarking
691 return nil, errors.New("not implemented, yet") 703 if err := generateMarkingPoints(
704 ctx, tx, feedback,
705 xyz, removed, epsg,
706 id,
707 ); err != nil {
708 return nil, err
709 }
692 } 710 }
693 711
694 // Store for potential later removal. 712 // Store for potential later removal.
695 if err = track(ctx, tx, importID, "waterway.sounding_results", id); err != nil { 713 if err = track(ctx, tx, importID, "waterway.sounding_results", id); err != nil {
696 return nil, err 714 return nil, err
906 heights = append(heights, h) 924 heights = append(heights, h)
907 } 925 }
908 return heights 926 return heights
909 } 927 }
910 928
929 func loadClassBreaks(
930 ctx context.Context,
931 tx *sql.Tx,
932 feedback Feedback,
933 minZ, maxZ float64,
934 ) mesh.ClassBreaks {
935
936 heights, err := mesh.LoadClassBreaks(
937 ctx, tx,
938 "morphology_classbreaks")
939
940 if err != nil {
941 feedback.Warn("Loading class breaks failed: %v", err)
942 feedback.Info("Using default class breaks")
943 heights = defaultClassBreaks(minZ, maxZ)
944 } else {
945 heights = heights.ExtrapolateClassBreaks(minZ, maxZ)
946 }
947
948 return heights.Dedup()
949 }
950
951 func generateMarkingPoints(
952 ctx context.Context,
953 tx *sql.Tx,
954 feedback Feedback,
955 xyz mesh.MultiPointZ,
956 removed map[int32]struct{},
957 epsg uint32,
958 id int64,
959 ) error {
960 log.Printf("debug: generateMarkingPoints")
961
962 min, max := mesh.MinMaxVertex(xyz.FilterRemoved(removed))
963
964 log.Printf("debug: min/max %.2f/%.2f\n", min.Z, max.Z)
965
966 heights := loadClassBreaks(ctx, tx, feedback, min.Z, max.Z)
967
968 classes := heights.Classify(xyz.FilterRemoved(removed))
969
970 log.Printf("debug: removed %d\n", len(removed))
971
972 // Should not happen ... Z values over the top.
973 if n := len(classes) - 1; n > 1 && len(classes[n]) > 0 {
974 // Place the over the top values to the class below.
975 classes[n-1] = append(classes[n-1], classes[n]...)
976 classes[n] = nil
977 classes = classes[:n]
978 }
979
980 stmt, err := tx.PrepareContext(ctx, insertMarkingPointsSQL)
981 if err != nil {
982 return err
983 }
984 defer stmt.Close()
985
986 for i, class := range classes {
987 // Ignore empty classes
988 if len(class) == 0 {
989 continue
990 }
991 log.Printf("debug: class %d: %d\n", i, len(class))
992 _, err := stmt.ExecContext(ctx, id, heights[i], epsg, class.AsWKB())
993 if err != nil {
994 return err
995 }
996 }
997
998 return nil
999 }
1000
911 func generateIsos( 1001 func generateIsos(
912 ctx context.Context, 1002 ctx context.Context,
913 tx *sql.Tx, 1003 tx *sql.Tx,
914 feedback Feedback, 1004 feedback Feedback,
915 tree *mesh.STRTree, 1005 tree *mesh.STRTree,
916 id int64, 1006 id int64,
917 ) error { 1007 ) error {
918 1008
919 heights, err := mesh.LoadClassBreaks( 1009 heights := loadClassBreaks(ctx, tx, feedback, tree.Min().Z, tree.Max().Z)
920 ctx, tx,
921 "morphology_classbreaks",
922 )
923
924 minZ, maxZ := tree.Min().Z, tree.Max().Z
925
926 if err != nil {
927 feedback.Warn("Loading class breaks failed: %v", err)
928 feedback.Info("Using default class breaks")
929 heights = defaultClassBreaks(minZ, maxZ)
930 } else {
931 heights = heights.ExtrapolateClassBreaks(minZ, maxZ)
932 }
933
934 /*
935 for i, v := range heights {
936 fmt.Printf("%d %.2f\n", i, v)
937 }
938 log.Printf("%.2f - %.2f\n", tree.Min.Z, tree.Max.Z)
939 */
940
941 heights = heights.Dedup()
942 1010
943 return generateIsoAreas(ctx, tx, feedback, tree, heights, id) 1011 return generateIsoAreas(ctx, tx, feedback, tree, heights, id)
944 } 1012 }
945 1013
946 func generateIsoAreas( 1014 func generateIsoAreas(