Mercurial > gemma

comparison pkg/controllers/cross.go @ 801:b6a1779ffb42

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Removed old cross section controller.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Thu, 27 Sep 2018 13:18:50 +0200
parents pkg/controllers/octreecross.go@0cc97135717c
children e10a8a682297
comparison
equal deleted inserted replaced
800:0689f4b7c99f 801:b6a1779ffb42
1 package controllers 1 package controllers
2 2
3 import ( 3 import (
4 "context"
4 "database/sql" 5 "database/sql"
5 "log" 6 "log"
6 "net/http" 7 "net/http"
7 "time" 8 "time"
8 9
9 "gemma.intevation.de/gemma/pkg/models" 10 "gemma.intevation.de/gemma/pkg/models"
11 "gemma.intevation.de/gemma/pkg/octree"
10 ) 12 )
11 13
12 const crossSQL = ` 14 const WGS84 = 4326
13 WITH line AS ( 15
14 SELECT ST_LineMerge(ST_Union(ST_3DIntersection( 16 func reproject(
15 ST_Translate( 17 rp *models.Reprojector,
16 ST_Extrude( 18 src models.GeoJSONLineCoordinates,
17 ST_GeomFromWKB($1, 4326), 19 ctx context.Context,
18 0, 0, 1000), 20 ) (models.GeoJSONLineCoordinates, error) {
19 0, 0, -500), 21
20 geom))) AS geom 22 dst := make(models.GeoJSONLineCoordinates, len(src))
21 FROM waterway.meshes m JOIN waterway.sounding_results sr ON m.sounding_result_id = sr.id 23 for i, s := range src {
22 WHERE ST_Intersects(geom, ST_GeomFromWKB($1, 4326)) AND 24 var err error
23 sr.bottleneck_id = $2 AND sr.date_info = $3 25 if dst[i].Lat, dst[i].Lon, err = rp.Reproject(
24 ) 26 s.Lat, s.Lon, ctx,
25 SELECT ST_AsBinary((ST_Dump(ST_Intersection(line.geom, sr.area::geometry))).geom) 27 ); err != nil {
26 FROM line, waterway.sounding_results sr 28 return nil, err
27 WHERE sr.bottleneck_id = $2 AND sr.date_info = $3 29 }
28 ` 30 }
31 return dst, nil
32 }
33
34 func projectBack(
35 rp *models.Reprojector,
36 lines octree.MultiLineStringZ,
37 ctx context.Context,
38 ) (models.GeoJSONMultiLineCoordinatesZ, error) {
39
40 out := make(models.GeoJSONMultiLineCoordinatesZ, len(lines))
41
42 for i, segment := range lines {
43
44 coords := make(models.GeoJSONLineCoordinatesZ, len(segment))
45 out[i] = coords
46 for j, v := range segment {
47 lat, lon, err := rp.Reproject(v.X, v.Y, ctx)
48 if err != nil {
49 return nil, err
50 }
51 coords[j] = models.GeoJSONCoordinateZ{Lat: lat, Lon: lon, Z: v.Z}
52 }
53 }
54 return out, nil
55 }
29 56
30 func crossSection( 57 func crossSection(
31 input interface{}, 58 input interface{},
32 req *http.Request, 59 req *http.Request,
33 db *sql.Conn, 60 conn *sql.Conn,
34 ) (jr JSONResult, err error) { 61 ) (jr JSONResult, err error) {
35 62
36 csi := input.(*models.CrossSectionInput) 63 csi := input.(*models.CrossSectionInput)
37 64
38 start := time.Now() 65 start := time.Now()
39 66
40 var rows *sql.Rows 67 tree, err := octree.Cache.Get(
41 if rows, err = db.QueryContext( 68 csi.Properties.Bottleneck, csi.Properties.Date.Time,
42 req.Context(), 69 conn, req.Context())
43 crossSQL, 70
44 csi.Geometry.Coordinates.AsWKB(), 71 log.Printf("loading octree took: %s\n", time.Since(start))
45 csi.Properties.Bottleneck, 72 if err != nil {
46 csi.Properties.Date.Time, 73 return
74 }
75
76 // The coordinate system of the octree is an UTM projection.
77 // The input coordinates are in WGS84.
78 // So we need to reproject them.
79
80 start = time.Now()
81
82 var rp *models.Reprojector
83 if rp, err = models.NewReprojector(
84 conn, req.Context(),
85 WGS84, tree.EPSG,
47 ); err != nil { 86 ); err != nil {
48 return 87 return
49 } 88 }
50 defer rows.Close() 89 defer rp.Close()
51 90
52 var segments models.GeoJSONMultiLineCoordinatesZ 91 coords, err := reproject(rp, csi.Geometry.Coordinates, req.Context())
53 92
54 for rows.Next() { 93 log.Printf("transforming input coords took: %s\n", time.Since(start))
55 var segment models.GeoJSONLineCoordinatesZ 94 if err != nil {
56 if err = rows.Scan(&segment); err != nil {
57 return
58 }
59 segments = append(segments, segment)
60 }
61
62 if err = rows.Err(); err != nil {
63 return 95 return
64 } 96 }
65 97
66 log.Printf("query took: %v\n", time.Since(start))
67
68 start = time.Now() 98 start = time.Now()
69 99
70 joined := segments.Join() 100 var segments octree.MultiLineStringZ
71 101
72 log.Printf("joining took: %v\n", time.Since(start)) 102 for i := 0; i < len(coords)-1; i++ {
73 log.Printf("segments before/after: %d %d\n", len(segments), len(joined)) 103 c1 := &coords[i]
104 c2 := &coords[i+1]
74 105
75 /* 106 verticalLine := octree.NewVerticalLine(c1.Lat, c1.Lon, c2.Lat, c2.Lon)
76 if err2 := dumpProfile( 107
77 csi.Geometry.Coordinates[0], 108 var line octree.MultiLineStringZ
78 csi.Geometry.Coordinates[len(csi.Geometry.Coordinates)-1], 109 tree.Vertical(c1.Lat, c1.Lon, c2.Lat, c2.Lon, func(t *octree.Triangle) {
79 joined, 110 if ls := verticalLine.Intersection(t); len(ls) > 0 {
80 ); err2 != nil { 111 line = append(line, ls)
81 log.Printf("error: %v\n", err2) 112 }
113 })
114
115 if len(line) > 0 {
116 log.Printf("line length: %d\n", len(line))
117 // They are all on the segment (c1.Lat, c1.Lon) - (c2.Lat, c2.Lon).
118 // Sort them by project them on this line.
119 joined := line.JoinOnLine(c1.Lat, c1.Lon, c2.Lat, c2.Lon)
120 log.Printf("joined length: %d\n", len(joined))
121 segments = append(segments, joined...)
82 } 122 }
83 */ 123
124 }
125 log.Printf("octree traversal took: %s\n", time.Since(start))
126
127 // The result have to be WGS84. So project the result back.
128 start = time.Now()
129
130 rp.FromEPSG, rp.ToEPSG = tree.EPSG, WGS84
131
132 var joined models.GeoJSONMultiLineCoordinatesZ
133 joined, err = projectBack(rp, segments, req.Context())
134
135 log.Printf("projecting back took: %s\n", time.Since(start))
136 if err != nil {
137 return
138 }
84 139
85 jr = JSONResult{ 140 jr = JSONResult{
86 Result: &models.CrossSectionOutput{ 141 Result: &models.CrossSectionOutput{
87 Type: "Feature", 142 Type: "Feature",
88 Geometry: models.CrossSectionOutputGeometry{ 143 Geometry: models.CrossSectionOutputGeometry{
93 }, 148 },
94 } 149 }
95 150
96 return 151 return
97 } 152 }
98
99 /*
100 func dumpProfile(
101 start models.GeoJSONCoordinate,
102 stop models.GeoJSONCoordinate,
103 segments models.GeoJSONMultiLineCoordinatesZ,
104 ) error {
105
106 w, err := os.Create("/tmp/cross.txt")
107 if err != nil {
108 return err
109 }
110 defer w.Close()
111
112 begin := models.GeoJSONCoordinateZ{
113 Lat: start.Lat,
114 Lon: start.Lon,
115 }
116 end := models.GeoJSONCoordinateZ{
117 Lat: stop.Lat,
118 Lon: stop.Lon,
119 }
120
121 last := begin
122
123 var total float64
124
125 minz := 10000.0
126
127 for _, line := range segments {
128 for _, coord := range line {
129 if coord.Z < minz {
130 minz = coord.Z
131 }
132 total += last.Distance(coord)
133 last = coord
134 }
135 }
136
137 log.Printf("length: %.3f minz: %f\n", total, minz)
138
139 var pos float64
140
141 for i, line := range segments {
142 for j, coord := range line {
143 if i == 0 && j == 0 {
144 if !begin.Equals(coord) {
145 pos = begin.Distance(coord)
146 fmt.Fprintf(w, "%.3f %f\n", 0.0, 200.0)
147 fmt.Fprintf(w, "%.3f %f\n", pos-0.0001, 200.0)
148 fmt.Fprintf(w, "%.3f %f\n", pos, coord.Z-minz)
149 continue
150 }
151 } else if j == 0 {
152 fmt.Fprintf(w, "%.3f %f\n", pos+0.0001, 200.0)
153 fmt.Fprintf(w, "%.3f %f\n", pos+last.Distance(coord)-0.0001, 200.0)
154 continue
155 }
156 pos += last.Distance(coord)
157 fmt.Fprintf(w, "%.3f %f\n", pos, coord.Z-minz)
158 last = coord
159 }
160 }
161
162 if !last.Equals(end) {
163 fmt.Fprintf(w, "%.3f %f\n", pos+0.0001, 200.0)
164 fmt.Fprintf(w, "%.3f %f\n", pos+last.Distance(end), 200.0)
165 }
166
167 return nil
168 }
169 */