Mercurial > gemma

comparison pkg/imports/sr.go @ 4564:6b107d4e6810 iso-areas

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add more debug output to trace problem why reprojected iso areas are not in expected bounds.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Thu, 03 Oct 2019 19:51:53 +0200
parents a4042ab02e05
children 2c49a293f275
comparison
equal deleted inserted replaced
4563:a4042ab02e05 4564:6b107d4e6810
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 "runtime" 32 "runtime"
204 ST_Transform( 205 ST_Transform(
205 ST_Multi( 206 ST_Multi(
206 ST_CollectionExtract( 207 ST_CollectionExtract(
207 ST_SimplifyPreserveTopology( 208 ST_SimplifyPreserveTopology(
208 ST_Multi(ST_Collectionextract( 209 ST_Multi(ST_Collectionextract(
209 ST_MakeValid(ST_GeomFromWKB($4, $3::integer)), 2)), 210 ST_MakeValid(ST_GeomFromWKB($4, $3::integer)), 3)),
210 $5 211 $5
211 ), 212 ),
212 3 213 3
213 ) 214 )
214 ), 215 ),
927 min, max := tree.Min, tree.Max 928 min, max := tree.Min, tree.Max
928 929
929 width := max.X - min.X 930 width := max.X - min.X
930 height := max.Y - min.Y 931 height := max.Y - min.Y
931 932
932 feedback.Info("width/height: %.2f / %.2f", width, height) 933 feedback.Info("Width/Height: %.2f / %.2f", width, height)
933 934
934 xcells := int(math.Ceil(width / isoCellSize)) 935 xcells := int(math.Ceil(width / isoCellSize))
935 ycells := int(math.Ceil(height / isoCellSize)) 936 ycells := int(math.Ceil(height / isoCellSize))
936 937
937 feedback.Info("raster size: (%d, %d)", xcells, ycells) 938 feedback.Info("Raster size: (%d, %d)", xcells, ycells)
938 939
939 // Add border for closing 940 // Add border for closing
940 raster := make([]float64, (xcells+2)*(ycells+2)) 941 raster := make([]float64, (xcells+2)*(ycells+2))
941 942
942 // prefill for no data 943 // prefill for no data
980 } 981 }
981 close(rows) 982 close(rows)
982 983
983 wg.Wait() 984 wg.Wait()
984 985
985 feedback.Info("Rasterizing took %v", time.Since(start)) 986 feedback.Info("Rasterizing took %v.", time.Since(start))
986 987
987 // TODO: Implement me! 988 feedback.Info("Trace iso areas.")
988 feedback.Info("Trace iso areas")
989 989
990 start = time.Now() 990 start = time.Now()
991 991
992 tracer := contourmap.FromFloat64s(xcells+2, ycells+2, raster) 992 tracer := contourmap.FromFloat64s(xcells+2, ycells+2, raster)
993 993
1002 doContours := func() { 1002 doContours := func() {
1003 defer wg.Done() 1003 defer wg.Done()
1004 for hIdx := range cnts { 1004 for hIdx := range cnts {
1005 c := tracer.Contours(heights[hIdx]) 1005 c := tracer.Contours(heights[hIdx])
1006 1006
1007 if len(c) == 0 {
1008 continue
1009 }
1010
1007 // We need to bring it back to the 1011 // We need to bring it back to the
1008 // none raster coordinate system. 1012 // none raster coordinate system.
1009 a := make(wkb.MultiPolygonGeom, len(c)) 1013 a := make(wkb.MultiPolygonGeom, len(c))
1014
1015 rXMin, rXMax := math.MaxFloat64, -math.MaxFloat64
1016 rYMin, rYMax := math.MaxFloat64, -math.MaxFloat64
1017
1010 for i, pl := range c { 1018 for i, pl := range c {
1011 shell := make(wkb.LinearRingGeom, len(pl)) 1019 shell := make(wkb.LinearRingGeom, len(pl))
1012 for j, pt := range pl { 1020 for j, pt := range pl {
1021 x := reprojX(pt.X)
1022 y := reprojY(pt.Y)
1023
1024 if x < rXMin {
1025 rXMin = x
1026 }
1027 if y < rYMin {
1028 rYMin = y
1029 }
1030 if x > rXMax {
1031 rXMax = x
1032 }
1033 if y > rYMax {
1034 rYMax = y
1035 }
1036
1013 shell[j] = wkb.PointGeom{ 1037 shell[j] = wkb.PointGeom{
1014 X: reprojX(pt.X), 1038 X: x,
1015 Y: reprojY(pt.Y), 1039 Y: y,
1016 } 1040 }
1017 } 1041 }
1018 a[i] = wkb.PolygonGeom{shell} 1042 a[i] = wkb.PolygonGeom{shell}
1019 } 1043 }
1044 log.Printf("(%.2f, %.2f) - (%.2f, %.2f)\n",
1045 rXMin, rYMin, rXMax, rYMax)
1046 log.Printf("%.2f / %.2f\n", rXMax-rXMin, rYMax-rYMin)
1047
1020 areas[hIdx] = a 1048 areas[hIdx] = a
1021 } 1049 }
1022 } 1050 }
1051
1052 log.Printf("[%.2f, %.2f] - [%.2f, %.2f]\n",
1053 min.X, min.Y, max.X, max.Y)
1023 1054
1024 for n := runtime.NumCPU(); n >= 1; n-- { 1055 for n := runtime.NumCPU(); n >= 1; n-- {
1025 wg.Add(1) 1056 wg.Add(1)
1026 go doContours() 1057 go doContours()
1027 } 1058 }
1031 } 1062 }
1032 close(cnts) 1063 close(cnts)
1033 1064
1034 wg.Wait() 1065 wg.Wait()
1035 1066
1036 feedback.Info("Tracing areas took %v", time.Since(start)) 1067 feedback.Info("Tracing areas took %v.", time.Since(start))
1037 1068
1038 // Raster and tracer are not needed any more. 1069 // Raster and tracer are not needed any more.
1039 raster, tracer = nil, nil 1070 raster, tracer = nil, nil
1040 1071
1041 return storeAreas( 1072 return storeAreas(
1050 areas []wkb.MultiPolygonGeom, 1081 areas []wkb.MultiPolygonGeom,
1051 epsg uint32, 1082 epsg uint32,
1052 heights []float64, 1083 heights []float64,
1053 id int64, 1084 id int64,
1054 ) error { 1085 ) error {
1055 feedback.Info("Storing iso areas") 1086 feedback.Info("Store iso areas.")
1056 total := time.Now() 1087 total := time.Now()
1057 defer func() { 1088 defer func() {
1058 feedback.Info("Storing iso areas took %v", 1089 feedback.Info("Storing iso areas took %v.",
1059 time.Since(total)) 1090 time.Since(total))
1060 }() 1091 }()
1092
1093 log.Printf("EPSG: %d\n", epsg)
1061 1094
1062 stmt, err := tx.PrepareContext(ctx, insertIsoAreasSQL) 1095 stmt, err := tx.PrepareContext(ctx, insertIsoAreasSQL)
1063 if err != nil { 1096 if err != nil {
1064 return err 1097 return err
1065 } 1098 }