Mercurial > gemma
diff cmd/octreediff/main.go @ 2479:c85b16db8a02 octree-diff
Calculate better triangulation and store it into database.
A bit fake by now. Needs extra:
CREATE TABLE diff_contour_lines (
height numeric NOT NULL,
lines geography(MultiLineString,4326) NOT NULL
);
GRANT ALL ON diff_contour_lines TO sophie;
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Thu, 28 Feb 2019 15:55:32 +0100 |
parents | efe332e985b9 |
children | 242104c338ff |
line wrap: on
line diff
--- a/cmd/octreediff/main.go Thu Feb 28 12:34:47 2019 +0100 +++ b/cmd/octreediff/main.go Thu Feb 28 15:55:32 2019 +0100 @@ -38,6 +38,8 @@ second = flag.String("second", "", "date of the second sounding result") ) +const contourTolerance = 0.1 + const ( loadOctreeSQL = ` SELECT sr.octree_index @@ -47,8 +49,56 @@ AND sr.octree_index IS NOT NULL` triangulateSQL = ` -SELECT ST_AsBinary(ST_DelaunayTriangles(ST_GeomFromWKB($1, $2::int), 0, 2)) +WITH joined AS ( + SELECT + sr.area AS area, + sr.date_info AS date_info + FROM waterway.sounding_results sr JOIN + waterway.bottlenecks bn ON sr.bottleneck_id = bn.id + WHERE bn.bottleneck_id = $1 +), +inter AS ( + SELECT + ST_Buffer( + ST_intersection( + (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $3::date), + (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $4::date) + ), + 0.001) AS area +), +triangles AS ( + SELECT t.geom AS geom, ST_MakePolygon(ST_ExteriorRing(t.geom)) AS poly FROM ( + SELECT (ST_Dump( + ST_DelaunayTriangles(ST_GeomFromWKB($5, $2::int), 0, 2))).geom + ) t +) +SELECT ST_AsBinary(ST_Collect(triangles.geom)) FROM triangles, inter +WHERE ST_Covers(inter.area, triangles.poly) ` + + //INSERT INTO redis_diff_countour_lines ( + //INSERT INTO diff_contour_lines ( + //INSERT INTO diff_contour_lines_extern ( + insertContourSQL = ` +INSERT INTO diff_contour_lines ( + height, + lines +) +SELECT + $1, + ST_Transform( + ST_Multi( + ST_CollectionExtract( + ST_SimplifyPreserveTopology( + ST_Multi(ST_Collectionextract( + ST_MakeValid(ST_GeomFromWKB($2, $3::integer)), 2)), + $4 + ), + 2 + ) + ), + 4326 + )` ) func check(err error) { @@ -311,8 +361,11 @@ if err := conn.QueryRowContext( ctx, triangulateSQL, + bottleneck, + first.EPSG, + firstDate, + secondDate, data, - first.EPSG, ).Scan(&tin); err != nil { return err } @@ -361,12 +414,31 @@ var dataSize int + stmt, err := conn.PrepareContext(ctx, insertContourSQL) + if err != nil { + return err + } + defer stmt.Close() + octree.DoContours(tree, heights, func(res *octree.ContourResult) { - // TODO: Store them. - log.Printf("%f: segments: %d\n", res.Height, len(res.Lines)) - dataSize += len(res.Lines.AsWKB2D()) + if err == nil && len(res.Lines) > 0 { + log.Printf("%f: lines: %d\n", res.Height, len(res.Lines)) + wkb := res.Lines.AsWKB2D() + dataSize += len(wkb) + _, err = stmt.ExecContext( + ctx, + res.Height, + wkb, + first.EPSG, + contourTolerance, + ) + } }) + if err != nil { + return err + } + now = time.Now() log.Printf("Number of iso lines: %d\n", len(heights)) log.Printf("Total WKB size: %.2fMB\n", float64(dataSize)/(1024*1024))