Mercurial > gemma

comparison pkg/imports/sr.go @ 4563:a4042ab02e05 iso-areas

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Store geometries in database.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Thu, 03 Oct 2019 10:05:01 +0200
parents 5cc4042cf07c
children 6b107d4e6810
comparison
equal deleted inserted replaced
4562:5cc4042cf07c 4563:a4042ab02e05
98 98
99 func (srJobCreator) Create() Job { return new(SoundingResult) } 99 func (srJobCreator) Create() Job { return new(SoundingResult) }
100 100
101 func (srJobCreator) Depends() [2][]string { 101 func (srJobCreator) Depends() [2][]string {
102 return [2][]string{ 102 return [2][]string{
103 {"sounding_results", "sounding_results_contour_lines"}, 103 {"sounding_results",
104 "sounding_results_contour_lines",
105 "sounding_results_iso_areas",
106 },
104 {"bottlenecks"}, 107 {"bottlenecks"},
105 } 108 }
106 } 109 }
107 110
108 func (srJobCreator) StageDone( 111 func (srJobCreator) StageDone(
204 ST_SimplifyPreserveTopology( 207 ST_SimplifyPreserveTopology(
205 ST_Multi(ST_Collectionextract( 208 ST_Multi(ST_Collectionextract(
206 ST_MakeValid(ST_GeomFromWKB($4, $3::integer)), 2)), 209 ST_MakeValid(ST_GeomFromWKB($4, $3::integer)), 2)),
207 $5 210 $5
208 ), 211 ),
209 2 212 3
210 ) 213 )
211 ), 214 ),
212 4326 215 4326
213 ) 216 )
214 FROM waterway.sounding_results sr 217 FROM waterway.sounding_results sr
988 991
989 tracer := contourmap.FromFloat64s(xcells+2, ycells+2, raster) 992 tracer := contourmap.FromFloat64s(xcells+2, ycells+2, raster)
990 993
991 areas := make([]wkb.MultiPolygonGeom, len(heights)) 994 areas := make([]wkb.MultiPolygonGeom, len(heights))
992 995
996 // TODO: Check if this correct!
993 reprojX := common.Linear(1, min.X, float64(xcells+1), max.X) 997 reprojX := common.Linear(1, min.X, float64(xcells+1), max.X)
994 reprojY := common.Linear(1, min.Y, float64(ycells+1), max.Y) 998 reprojY := common.Linear(1, min.Y, float64(ycells+1), max.Y)
995 999
996 cnts := make(chan int) 1000 cnts := make(chan int)
997 1001
998 doContours := func() { 1002 doContours := func() {
999 defer wg.Done() 1003 defer wg.Done()
1000 for hIdx := range cnts { 1004 for hIdx := range cnts {
1001 c := tracer.Contours(heights[hIdx]) 1005 c := tracer.Contours(heights[hIdx])
1006
1007 // We need to bring it back to the
1008 // none raster coordinate system.
1002 a := make(wkb.MultiPolygonGeom, len(c)) 1009 a := make(wkb.MultiPolygonGeom, len(c))
1003
1004 for i, pl := range c { 1010 for i, pl := range c {
1005
1006 shell := make(wkb.LinearRingGeom, len(pl)) 1011 shell := make(wkb.LinearRingGeom, len(pl))
1007 for j, pt := range pl { 1012 for j, pt := range pl {
1008 x := reprojX(pt.X) 1013 shell[j] = wkb.PointGeom{
1009 y := reprojY(pt.Y) 1014 X: reprojX(pt.X),
1010 shell[j] = wkb.PointGeom{X: x, Y: y} 1015 Y: reprojY(pt.Y),
1016 }
1011 } 1017 }
1012
1013 a[i] = wkb.PolygonGeom{shell} 1018 a[i] = wkb.PolygonGeom{shell}
1014 } 1019 }
1015 areas[hIdx] = a 1020 areas[hIdx] = a
1016 } 1021 }
1017 } 1022 }
1031 feedback.Info("Tracing areas took %v", time.Since(start)) 1036 feedback.Info("Tracing areas took %v", time.Since(start))
1032 1037
1033 // Raster and tracer are not needed any more. 1038 // Raster and tracer are not needed any more.
1034 raster, tracer = nil, nil 1039 raster, tracer = nil, nil
1035 1040
1036 return storeAreas(ctx, tx, feedback, areas, heights) 1041 return storeAreas(
1042 ctx, tx, feedback,
1043 areas, tree.EPSG, heights, id)
1037 } 1044 }
1038 1045
1039 func storeAreas( 1046 func storeAreas(
1040 ctx context.Context, 1047 ctx context.Context,
1041 tx *sql.Tx, 1048 tx *sql.Tx,
1042 feedback Feedback, 1049 feedback Feedback,
1043 areas []wkb.MultiPolygonGeom, 1050 areas []wkb.MultiPolygonGeom,
1051 epsg uint32,
1044 heights []float64, 1052 heights []float64,
1053 id int64,
1045 ) error { 1054 ) error {
1046 feedback.Info("Storing iso areas") 1055 feedback.Info("Storing iso areas")
1047 total := time.Now() 1056 total := time.Now()
1048 defer func() { 1057 defer func() {
1049 feedback.Info("Storing iso areas took %v", 1058 feedback.Info("Storing iso areas took %v",
1054 if err != nil { 1063 if err != nil {
1055 return err 1064 return err
1056 } 1065 }
1057 defer stmt.Close() 1066 defer stmt.Close()
1058 1067
1059 // TODO: Implement me! 1068 for i, a := range areas {
1069 if len(a) == 0 {
1070 continue
1071 }
1072 if _, err := stmt.ExecContext(
1073 ctx,
1074 id, heights[i], epsg,
1075 a.AsWKB(),
1076 contourTolerance,
1077 ); err != nil {
1078 return err
1079 }
1080 }
1060 1081
1061 return nil 1082 return nil
1062 } 1083 }
1063 1084
1064 func generateContours( 1085 func generateContours(