# HG changeset patch # User Sascha L. Teichmann # Date 1551891118 -3600 # Node ID 45d51a49f191606504e43e83292c4f23f5e642e7 # Parent 113912129481f5598972938d7cf565060cba87db SR import: Use own triangulation and clipping when importing sounding results. diff -r 113912129481 -r 45d51a49f191 pkg/imports/sr.go --- 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(), diff -r 113912129481 -r 45d51a49f191 pkg/octree/contours.go --- a/pkg/octree/contours.go Wed Mar 06 16:26:45 2019 +0100 +++ b/pkg/octree/contours.go Wed Mar 06 17:51:58 2019 +0100 @@ -15,7 +15,6 @@ package octree import ( - "log" "runtime" "sync" ) @@ -71,8 +70,6 @@ return } - log.Printf("num heights: %d\n", len(heights)) - contours := make([]*ContourResult, len(heights)) for i, h := range heights { diff -r 113912129481 -r 45d51a49f191 pkg/octree/tin.go --- a/pkg/octree/tin.go Wed Mar 06 16:26:45 2019 +0100 +++ b/pkg/octree/tin.go Wed Mar 06 17:51:58 2019 +0100 @@ -15,8 +15,6 @@ import ( "bytes" - "context" - "database/sql" "encoding/binary" "errors" "fmt" @@ -190,47 +188,6 @@ return nil } -const ( - tinSQLPrefix = `WITH trans AS ( - SELECT - ST_Buffer(ST_Transform(area::geometry, $1::int), 0.001) AS area, - ST_Transform(point_cloud::geometry, $1::int) AS point_cloud - FROM waterway.sounding_results -` - tinSQLSuffix = ` -), -triangles AS ( - SELECT t.geom AS geom, ST_MakePolygon(ST_ExteriorRing(t.geom)) AS poly FROM ( - SELECT (ST_Dump( - ST_DelaunayTriangles(point_cloud, 0, 2))).geom - FROM trans) t -) -SELECT ST_AsBinary(ST_Collect(triangles.geom)) FROM triangles, trans -WHERE ST_Covers(trans.area, triangles.poly)` - - loadTinByIDSQL = tinSQLPrefix + `WHERE id = $2` + tinSQLSuffix -) - -// GenerateTinByID generated a TIN by triangulating a point cloud -// from the database. -func GenerateTinByID( - ctx context.Context, - conn *sql.Conn, - id int64, - epsg uint32, -) (*Tin, error) { - var tin Tin - err := conn.QueryRowContext(ctx, loadTinByIDSQL, epsg, id).Scan(&tin) - switch { - case err == sql.ErrNoRows: - return nil, nil - case err != nil: - return nil, err - } - tin.EPSG = epsg - return &tin, nil -} - // Scan implements the sql.Scanner interface. func (t *Tin) Scan(raw interface{}) error { if raw == nil { diff -r 113912129481 -r 45d51a49f191 pkg/octree/vertex.go --- a/pkg/octree/vertex.go Wed Mar 06 16:26:45 2019 +0100 +++ b/pkg/octree/vertex.go Wed Mar 06 17:51:58 2019 +0100 @@ -16,6 +16,7 @@ import ( "bytes" "encoding/binary" + "fmt" "io" "log" "math" @@ -1059,3 +1060,84 @@ return buf.Bytes() } + +func (mpz *MultiPointZ) FromWKB(data []byte) error { + + r := bytes.NewReader(data) + + endian, err := r.ReadByte() + + var order binary.ByteOrder + + switch { + case err != nil: + return err + case endian == wkb.NDR: + order = binary.LittleEndian + case endian == wkb.XDR: + order = binary.BigEndian + default: + return fmt.Errorf("unknown byte order %x", endian) + } + + var geomType uint32 + err = binary.Read(r, order, &geomType) + + switch { + case err != nil: + return err + case geomType != wkb.MultiPointZ: + return fmt.Errorf("unknown geometry type %x", geomType) + } + + var numPoints uint32 + if err := binary.Read(r, order, &numPoints); err != nil { + return err + } + + points := make(MultiPointZ, numPoints) + + for i := range points { + endian, err = r.ReadByte() + + switch { + case err != nil: + return err + case endian == wkb.NDR: + order = binary.LittleEndian + case endian == wkb.XDR: + order = binary.BigEndian + default: + return fmt.Errorf("unknown byte order %x", endian) + } + + err = binary.Read(r, order, &geomType) + + switch { + case err != nil: + return err + case geomType != wkb.PointZ: + return fmt.Errorf("unknown geometry type %x", geomType) + } + + var x, y, z uint64 + if err = binary.Read(r, order, &x); err != nil { + return err + } + if err = binary.Read(r, order, &y); err != nil { + return err + } + if err = binary.Read(r, order, &z); err != nil { + return err + } + points[i] = Vertex{ + X: math.Float64frombits(x), + Y: math.Float64frombits(y), + Z: math.Float64frombits(z), + } + } + + *mpz = points + + return nil +} diff -r 113912129481 -r 45d51a49f191 schema/gemma.sql --- a/schema/gemma.sql Wed Mar 06 16:26:45 2019 +0100 +++ b/schema/gemma.sql Wed Mar 06 17:51:58 2019 +0100 @@ -505,8 +505,6 @@ surtyp varchar REFERENCES survey_types, coverage varchar REFERENCES coverage_types, depth_reference varchar(4) NOT NULL, -- REFERENCES depth_references, - point_cloud geography(MULTIPOINTZ, 4326) NOT NULL - CHECK(ST_IsSimple(CAST(point_cloud AS geometry))), octree_checksum varchar, octree_index bytea, staging_done boolean NOT NULL DEFAULT false