Mercurial > gemma
changeset 3600:6590733e2ebc
SR import: If in single beam mode and no boundary is given eliminate triangles which have at least an edge with is longer than 98% of the overall edge lengths.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Tue, 04 Jun 2019 18:37:15 +0200 |
parents | 73add36de413 |
children | 98f14d97611e |
files | pkg/imports/sr.go |
diffstat | 1 files changed, 58 insertions(+), 5 deletions(-) [+] |
line wrap: on
line diff
--- a/pkg/imports/sr.go Tue Jun 04 18:11:23 2019 +0200 +++ b/pkg/imports/sr.go Tue Jun 04 18:37:15 2019 +0200 @@ -28,11 +28,13 @@ "os" "path" "path/filepath" + "sort" "strconv" "strings" "time" shp "github.com/jonas-p/go-shp" + "gonum.org/v1/gonum/stat" "gemma.intevation.de/gemma/pkg/common" "gemma.intevation.de/gemma/pkg/misc" @@ -176,11 +178,12 @@ ` reprojectPointsSingleBeamSQL = ` -SELECT ST_AsBinary( - ST_Transform( - ST_GeomFromWKB($1, $2::integer), - best_utm(CAST(ST_GeomFromWKB($1, $2::integer) AS geography))), - best_utm(CAST(ST_GeomFromWKB($1, $2::integer) AS geography)) +SELECT + ST_AsBinary( + ST_Transform( + ST_GeomFromWKB($1, $2::integer), + best_utm(CAST(ST_Transform(ST_GeomFromWKB($1, $2::integer), 4326) AS geography)))), + best_utm(CAST(ST_Transform(ST_GeomFromWKB($1, $2::integer), 4326) AS geography)) ` ) @@ -205,6 +208,7 @@ defer z.Close() feedback.Info("Looking for 'meta.json'") + mf := common.FindInZIP(z, "meta.json") if mf == nil && !sr.completeOverride() { return nil, errors.New("Cannot find 'meta.json'") @@ -342,8 +346,57 @@ return nil, err } + if err := xyz.FromWKB(reproj); err != nil { + return nil, err + } + feedback.Info("Reprojecting points to EPSG %d took %v.", epsg, time.Since(start)) + feedback.Info("Number of reprojected points: %d", len(xyz)) + + start = time.Now() + tri, err := octree.Triangulate(xyz) + if err != nil { + return nil, err + } + feedback.Info("Triangulation took %v.", time.Since(start)) + feedback.Info("Number triangles: %d.", len(tri.Triangles)/3) + + if boundary == nil { + feedback.Info("No boundary given.") + feedback.Info("Eliminate triangles with long edges.") + + lens := make([]float64, len(tri.Triangles)/3) + for i := 0; i < len(tri.Triangles); i += 3 { + vs := tri.Triangles[i : i+3] + var max float64 + for j := 0; j < 3; j++ { + k := (j + 1) % 3 + if l := tri.Points[vs[j]].SquaredDistance2D(tri.Points[vs[k]]); l > max { + max = l + } + } + lens[i/3] = max + } + qs := make([]float64, len(lens)) + copy(qs, lens) + sort.Float64s(qs) + q := stat.Quantile(0.98, stat.LinInterp, qs, nil) + feedback.Info("98%% quantile of the edge lengths: %.2f\n", math.Sqrt(q)) + + ts := tri.Triangles + + for i := len(lens) - 1; i >= 0; i-- { + if lens[i] >= q { + if i != len(lens)-1 { + copy(ts[i*3:], ts[(i+1)*3:]) + } + ts = ts[:len(ts)-3] + } + } + feedback.Info("Remaining triangles: %d (%.2f%%)\n", + len(ts)/3, float64(len(ts))/float64(len(tri.Triangles))*100.0) + } // TODO: Implement me! return nil, errors.New("Not implemented, yet!")