Mercurial > gemma
diff pkg/imports/sr.go @ 2529:45d51a49f191
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 |
line wrap: on
line diff
--- a/pkg/imports/sr.go Wed Mar 06 16:26:45 2019 +0100 +++ b/pkg/imports/sr.go Wed Mar 06 17:51:58 2019 +0100 @@ -104,36 +104,33 @@ WHERE import_id = $1 AND relation = 'waterway.sounding_results'::regclass)` - insertPointsSQL = ` + insertHullSQL = ` INSERT INTO waterway.sounding_results ( bottleneck_id, date_info, depth_reference, - point_cloud, area -) VALUES ( +) SELECT (SELECT id from waterway.bottlenecks where objnam = $1), $2::date, $3, - ST_Transform(ST_GeomFromWKB($4, $6::integer), 4326)::geography, - (SELECT CASE WHEN $5::bytea IS NULL THEN - ST_Transform( - ST_ConcaveHull( - ST_Force2D(ST_GeomFromWKB($4, $6::integer)), 0.7), 4326)::geography - ELSE - ST_Transform(ST_GeomFromWKB($5, $6::integer), 4326)::geography - END) -) + (SELECT + CASE WHEN $5::bytea IS NULL THEN + ST_Transform(ST_ConcaveHull(ST_Force2D(ST_GeomFromWKB($4, $6::integer)), 0.7), 4326)::geography + ELSE + ST_Transform(ST_GeomFromWKB($5, $6::integer), 4326)::geography + END) RETURNING id, - ST_X(ST_Centroid(point_cloud::geometry)), - ST_Y(ST_Centroid(point_cloud::geometry)), - CASE WHEN ST_Y(ST_Centroid(point_cloud::geometry)) > 0 THEN - 32600 - ELSE - 32700 - END + floor((ST_X(ST_Centroid(point_cloud::geometry))+180)/6)::int + 1` + ST_X(ST_Centroid(area::geometry)), + ST_Y(ST_Centroid(area::geometry)), + best_utm(area::geometry), + ST_AsBinary(ST_Transform(area::geometry, best_utm(area::geometry))) +` + reprojectPointsSQL = ` +SELECT ST_AsBinary(ST_Transform(ST_GeomFromWKB($1, $2::integer), $3::integer)) +` insertOctreeSQL = ` UPDATE waterway.sounding_results SET octree_checksum = $2, octree_index = $3 @@ -174,6 +171,8 @@ feedback Feedback, ) (interface{}, error) { + begin := time.Now() + z, err := zip.OpenReader(filepath.Join(sr.Dir, "sr.zip")) if err != nil { return nil, err @@ -234,63 +233,116 @@ ) start := time.Now() - err = tx.QueryRow(insertPointsSQL, + var hull []byte + + xyzWKB := xyz.AsWKB() + + err = tx.QueryRow(insertHullSQL, m.Bottleneck, m.Date.Time, m.DepthReference, - xyz.AsWKB(), + xyzWKB, polygon.asWKB(), m.EPSG, - ).Scan(&id, &lat, &lon, &epsg) + ).Scan( + &id, + &lat, + &lon, + &epsg, + &hull, + ) xyz, polygon = nil, nil // not need from now on. - feedback.Info("storing points took %s", time.Since(start)) + feedback.Info("Calculating hull took %s.", time.Since(start)) + if err != nil { + return nil, err + } + feedback.Info("Best suited UTM EPSG: %d", epsg) + + start = time.Now() + + var clippingPolygon octree.Polygon + if err := clippingPolygon.FromWKB(hull); err != nil { + return nil, err + } + clippingPolygon.Indexify() + + feedback.Info("Building clipping polygon took %v.", time.Since(start)) + + start = time.Now() + + var reproj []byte + + err = tx.QueryRow(reprojectPointsSQL, + xyzWKB, + m.EPSG, + epsg, + ).Scan(&reproj) + if err != nil { return nil, err } - feedback.Info("Best suited UTM EPSG: %d", epsg) + if err := xyz.FromWKB(reproj); err != nil { + return nil, err + } - feedback.Info("Triangulate...") + feedback.Info("Reprojecting points took %v.", time.Since(start)) + feedback.Info("Number of reprojected points: %d", len(xyz)) + start = time.Now() - tin, err := octree.GenerateTinByID(ctx, conn, id, epsg) - feedback.Info("triangulation took %s", time.Since(start)) + + tri, err := octree.Triangulate(xyz) if err != nil { return nil, err } - if tin == nil { - return nil, errors.New("cannot load TIN from database") - } + feedback.Info("Triangulation took %v.", time.Since(start)) + + start = time.Now() + + tin := tri.Tin() + + var str octree.STRTree + str.Build(tin) - feedback.Info("Building octree...") + feedback.Info("Building STR tree took %v", time.Since(start)) + + start = time.Now() + + removed := str.Clip(&clippingPolygon) + feedback.Info("Clipping STR tree took %v.", time.Since(start)) + feedback.Info("Number of triangles to clip %d.", len(removed)) + + start = time.Now() + builder := octree.NewBuilder(tin) - start = time.Now() - builder.Build(nil) + builder.Build(removed) octreeIndex, err := builder.Bytes() - tin = nil // not needed from now on - feedback.Info("building octree took %s", time.Since(start)) + feedback.Info("Building octree took %v.", + time.Since(start)) if err != nil { return nil, err } - feedback.Info("Store octree...") start = time.Now() h := sha1.New() h.Write(octreeIndex) checksum := hex.EncodeToString(h.Sum(nil)) _, err = tx.Exec(insertOctreeSQL, id, checksum, octreeIndex) - feedback.Info("storing octree index took %s", time.Since(start)) + feedback.Info("Storing octree index took %s.", + time.Since(start)) if err != nil { return nil, err } tree := builder.Tree() + tree.EPSG = epsg builder = nil // not needed from now on - feedback.Info("Generate contour lines...") start = time.Now() err = generateContours(tree, tx, id) - feedback.Info("generating and storing contour lines took %s", time.Since(start)) + feedback.Info("Generating and storing contour lines took %s.", + time.Since(start)) if err != nil { return nil, err } @@ -300,10 +352,17 @@ return nil, err } - if err = tx.Commit(); err == nil { - feedback.Info("Storing sounding result was successful.") + if err = tx.Commit(); err != nil { + feedback.Error( + "Storing sounding result failed after %v.", + time.Since(begin)) + return nil, err } + feedback.Info( + "Storing sounding result was successful after %v.", + time.Since(begin)) + summary := struct { Bottleneck string `json:"bottleneck"` Date models.Date `json:"date"` @@ -473,7 +532,7 @@ } octree.DoContours(tree, heights, func(res *octree.ContourResult) { - if err == nil { + if err == nil && len(res.Lines) > 0 { _, err = stmt.Exec( id, res.Height, tree.EPSG, res.Lines.AsWKB2D(),