Mercurial > gemma

comparison pkg/imports/sr.go @ 2529:45d51a49f191

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
SR import: Use own triangulation and clipping when importing sounding results.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Wed, 06 Mar 2019 17:51:58 +0100
parents 1ec4c5633eb6
children 3f61b84ae7a6
comparison
equal deleted inserted replaced
2528:113912129481 2529:45d51a49f191
102 WHERE id = ( 102 WHERE id = (
103 SELECT key from import.track_imports 103 SELECT key from import.track_imports
104 WHERE import_id = $1 AND 104 WHERE import_id = $1 AND
105 relation = 'waterway.sounding_results'::regclass)` 105 relation = 'waterway.sounding_results'::regclass)`
106 106
107 insertPointsSQL = ` 107 insertHullSQL = `
108 INSERT INTO waterway.sounding_results ( 108 INSERT INTO waterway.sounding_results (
109 bottleneck_id, 109 bottleneck_id,
110 date_info, 110 date_info,
111 depth_reference, 111 depth_reference,
112 point_cloud,
113 area 112 area
114 ) VALUES ( 113 ) SELECT
115 (SELECT id from waterway.bottlenecks where objnam = $1), 114 (SELECT id from waterway.bottlenecks where objnam = $1),
116 $2::date, 115 $2::date,
117 $3, 116 $3,
118 ST_Transform(ST_GeomFromWKB($4, $6::integer), 4326)::geography, 117 (SELECT
119 (SELECT CASE WHEN $5::bytea IS NULL THEN 118 CASE WHEN $5::bytea IS NULL THEN
120 ST_Transform( 119 ST_Transform(ST_ConcaveHull(ST_Force2D(ST_GeomFromWKB($4, $6::integer)), 0.7), 4326)::geography
121 ST_ConcaveHull( 120 ELSE
122 ST_Force2D(ST_GeomFromWKB($4, $6::integer)), 0.7), 4326)::geography 121 ST_Transform(ST_GeomFromWKB($5, $6::integer), 4326)::geography
123 ELSE 122 END)
124 ST_Transform(ST_GeomFromWKB($5, $6::integer), 4326)::geography
125 END)
126 )
127 RETURNING 123 RETURNING
128 id, 124 id,
129 ST_X(ST_Centroid(point_cloud::geometry)), 125 ST_X(ST_Centroid(area::geometry)),
130 ST_Y(ST_Centroid(point_cloud::geometry)), 126 ST_Y(ST_Centroid(area::geometry)),
131 CASE WHEN ST_Y(ST_Centroid(point_cloud::geometry)) > 0 THEN 127 best_utm(area::geometry),
132 32600 128 ST_AsBinary(ST_Transform(area::geometry, best_utm(area::geometry)))
133 ELSE 129 `
134 32700 130
135 END + floor((ST_X(ST_Centroid(point_cloud::geometry))+180)/6)::int + 1` 131 reprojectPointsSQL = `
136 132 SELECT ST_AsBinary(ST_Transform(ST_GeomFromWKB($1, $2::integer), $3::integer))
133 `
137 insertOctreeSQL = ` 134 insertOctreeSQL = `
138 UPDATE waterway.sounding_results SET 135 UPDATE waterway.sounding_results SET
139 octree_checksum = $2, octree_index = $3 136 octree_checksum = $2, octree_index = $3
140 WHERE id = $1` 137 WHERE id = $1`
141 138
172 importID int64, 169 importID int64,
173 conn *sql.Conn, 170 conn *sql.Conn,
174 feedback Feedback, 171 feedback Feedback,
175 ) (interface{}, error) { 172 ) (interface{}, error) {
176 173
174 begin := time.Now()
175
177 z, err := zip.OpenReader(filepath.Join(sr.Dir, "sr.zip")) 176 z, err := zip.OpenReader(filepath.Join(sr.Dir, "sr.zip"))
178 if err != nil { 177 if err != nil {
179 return nil, err 178 return nil, err
180 } 179 }
181 defer z.Close() 180 defer z.Close()
232 epsg uint32 231 epsg uint32
233 lat, lon float64 232 lat, lon float64
234 ) 233 )
235 start := time.Now() 234 start := time.Now()
236 235
237 err = tx.QueryRow(insertPointsSQL, 236 var hull []byte
237
238 xyzWKB := xyz.AsWKB()
239
240 err = tx.QueryRow(insertHullSQL,
238 m.Bottleneck, 241 m.Bottleneck,
239 m.Date.Time, 242 m.Date.Time,
240 m.DepthReference, 243 m.DepthReference,
241 xyz.AsWKB(), 244 xyzWKB,
242 polygon.asWKB(), 245 polygon.asWKB(),
243 m.EPSG, 246 m.EPSG,
244 ).Scan(&id, &lat, &lon, &epsg) 247 ).Scan(
248 &id,
249 &lat,
250 &lon,
251 &epsg,
252 &hull,
253 )
245 xyz, polygon = nil, nil // not need from now on. 254 xyz, polygon = nil, nil // not need from now on.
246 feedback.Info("storing points took %s", time.Since(start)) 255 feedback.Info("Calculating hull took %s.", time.Since(start))
247 if err != nil { 256 if err != nil {
248 return nil, err 257 return nil, err
249 } 258 }
250
251 feedback.Info("Best suited UTM EPSG: %d", epsg) 259 feedback.Info("Best suited UTM EPSG: %d", epsg)
252 260
253 feedback.Info("Triangulate...") 261 start = time.Now()
254 start = time.Now() 262
255 tin, err := octree.GenerateTinByID(ctx, conn, id, epsg) 263 var clippingPolygon octree.Polygon
256 feedback.Info("triangulation took %s", time.Since(start)) 264 if err := clippingPolygon.FromWKB(hull); err != nil {
257 if err != nil { 265 return nil, err
258 return nil, err 266 }
259 } 267 clippingPolygon.Indexify()
260 268
261 if tin == nil { 269 feedback.Info("Building clipping polygon took %v.", time.Since(start))
262 return nil, errors.New("cannot load TIN from database") 270
263 } 271 start = time.Now()
264 272
265 feedback.Info("Building octree...") 273 var reproj []byte
274
275 err = tx.QueryRow(reprojectPointsSQL,
276 xyzWKB,
277 m.EPSG,
278 epsg,
279 ).Scan(&reproj)
280
281 if err != nil {
282 return nil, err
283 }
284
285 if err := xyz.FromWKB(reproj); err != nil {
286 return nil, err
287 }
288
289 feedback.Info("Reprojecting points took %v.", time.Since(start))
290 feedback.Info("Number of reprojected points: %d", len(xyz))
291
292 start = time.Now()
293
294 tri, err := octree.Triangulate(xyz)
295 if err != nil {
296 return nil, err
297 }
298
299 feedback.Info("Triangulation took %v.", time.Since(start))
300
301 start = time.Now()
302
303 tin := tri.Tin()
304
305 var str octree.STRTree
306 str.Build(tin)
307
308 feedback.Info("Building STR tree took %v", time.Since(start))
309
310 start = time.Now()
311
312 removed := str.Clip(&clippingPolygon)
313 feedback.Info("Clipping STR tree took %v.", time.Since(start))
314 feedback.Info("Number of triangles to clip %d.", len(removed))
315
316 start = time.Now()
317
266 builder := octree.NewBuilder(tin) 318 builder := octree.NewBuilder(tin)
267 start = time.Now() 319 builder.Build(removed)
268 builder.Build(nil)
269 octreeIndex, err := builder.Bytes() 320 octreeIndex, err := builder.Bytes()
270 tin = nil // not needed from now on 321 feedback.Info("Building octree took %v.",
271 feedback.Info("building octree took %s", time.Since(start)) 322 time.Since(start))
272 if err != nil { 323 if err != nil {
273 return nil, err 324 return nil, err
274 } 325 }
275 326
276 feedback.Info("Store octree...")
277 start = time.Now() 327 start = time.Now()
278 h := sha1.New() 328 h := sha1.New()
279 h.Write(octreeIndex) 329 h.Write(octreeIndex)
280 checksum := hex.EncodeToString(h.Sum(nil)) 330 checksum := hex.EncodeToString(h.Sum(nil))
281 _, err = tx.Exec(insertOctreeSQL, id, checksum, octreeIndex) 331 _, err = tx.Exec(insertOctreeSQL, id, checksum, octreeIndex)
282 feedback.Info("storing octree index took %s", time.Since(start)) 332 feedback.Info("Storing octree index took %s.",
333 time.Since(start))
283 if err != nil { 334 if err != nil {
284 return nil, err 335 return nil, err
285 } 336 }
286 337
287 tree := builder.Tree() 338 tree := builder.Tree()
339 tree.EPSG = epsg
288 builder = nil // not needed from now on 340 builder = nil // not needed from now on
289 341
290 feedback.Info("Generate contour lines...")
291 start = time.Now() 342 start = time.Now()
292 err = generateContours(tree, tx, id) 343 err = generateContours(tree, tx, id)
293 feedback.Info("generating and storing contour lines took %s", time.Since(start)) 344 feedback.Info("Generating and storing contour lines took %s.",
345 time.Since(start))
294 if err != nil { 346 if err != nil {
295 return nil, err 347 return nil, err
296 } 348 }
297 349
298 // Store for potential later removal. 350 // Store for potential later removal.
299 if err = track(ctx, tx, importID, "waterway.sounding_results", id); err != nil { 351 if err = track(ctx, tx, importID, "waterway.sounding_results", id); err != nil {
300 return nil, err 352 return nil, err
301 } 353 }
302 354
303 if err = tx.Commit(); err == nil { 355 if err = tx.Commit(); err != nil {
304 feedback.Info("Storing sounding result was successful.") 356 feedback.Error(
305 } 357 "Storing sounding result failed after %v.",
358 time.Since(begin))
359 return nil, err
360 }
361
362 feedback.Info(
363 "Storing sounding result was successful after %v.",
364 time.Since(begin))
306 365
307 summary := struct { 366 summary := struct {
308 Bottleneck string `json:"bottleneck"` 367 Bottleneck string `json:"bottleneck"`
309 Date models.Date `json:"date"` 368 Date models.Date `json:"date"`
310 Lat float64 `json:"lat"` 369 Lat float64 `json:"lat"`
471 for ; h <= tree.Max.Z; h += contourStepWidth { 530 for ; h <= tree.Max.Z; h += contourStepWidth {
472 heights = append(heights, h) 531 heights = append(heights, h)
473 } 532 }
474 533
475 octree.DoContours(tree, heights, func(res *octree.ContourResult) { 534 octree.DoContours(tree, heights, func(res *octree.ContourResult) {
476 if err == nil { 535 if err == nil && len(res.Lines) > 0 {
477 _, err = stmt.Exec( 536 _, err = stmt.Exec(
478 id, res.Height, tree.EPSG, 537 id, res.Height, tree.EPSG,
479 res.Lines.AsWKB2D(), 538 res.Lines.AsWKB2D(),
480 contourTolerance) 539 contourTolerance)
481 } 540 }