# HG changeset patch # User Sascha L. Teichmann # Date 1538047130 -7200 # Node ID b6a1779ffb42df8599472f1c9285eb0c376488f4 # Parent 0689f4b7c99f0c1112bea4c3d85d63751cea2437 Removed old cross section controller. diff -r 0689f4b7c99f -r b6a1779ffb42 pkg/controllers/cross.go --- a/pkg/controllers/cross.go Thu Sep 27 12:49:45 2018 +0200 +++ b/pkg/controllers/cross.go Thu Sep 27 13:18:50 2018 +0200 @@ -1,86 +1,141 @@ package controllers import ( + "context" "database/sql" "log" "net/http" "time" "gemma.intevation.de/gemma/pkg/models" + "gemma.intevation.de/gemma/pkg/octree" ) -const crossSQL = ` -WITH line AS ( -SELECT ST_LineMerge(ST_Union(ST_3DIntersection( - ST_Translate( - ST_Extrude( - ST_GeomFromWKB($1, 4326), - 0, 0, 1000), - 0, 0, -500), - geom))) AS geom -FROM waterway.meshes m JOIN waterway.sounding_results sr ON m.sounding_result_id = sr.id -WHERE ST_Intersects(geom, ST_GeomFromWKB($1, 4326)) AND - sr.bottleneck_id = $2 AND sr.date_info = $3 -) -SELECT ST_AsBinary((ST_Dump(ST_Intersection(line.geom, sr.area::geometry))).geom) - FROM line, waterway.sounding_results sr - WHERE sr.bottleneck_id = $2 AND sr.date_info = $3 -` +const WGS84 = 4326 + +func reproject( + rp *models.Reprojector, + src models.GeoJSONLineCoordinates, + ctx context.Context, +) (models.GeoJSONLineCoordinates, error) { + + dst := make(models.GeoJSONLineCoordinates, len(src)) + for i, s := range src { + var err error + if dst[i].Lat, dst[i].Lon, err = rp.Reproject( + s.Lat, s.Lon, ctx, + ); err != nil { + return nil, err + } + } + return dst, nil +} + +func projectBack( + rp *models.Reprojector, + lines octree.MultiLineStringZ, + ctx context.Context, +) (models.GeoJSONMultiLineCoordinatesZ, error) { + + out := make(models.GeoJSONMultiLineCoordinatesZ, len(lines)) + + for i, segment := range lines { + + coords := make(models.GeoJSONLineCoordinatesZ, len(segment)) + out[i] = coords + for j, v := range segment { + lat, lon, err := rp.Reproject(v.X, v.Y, ctx) + if err != nil { + return nil, err + } + coords[j] = models.GeoJSONCoordinateZ{Lat: lat, Lon: lon, Z: v.Z} + } + } + return out, nil +} func crossSection( input interface{}, req *http.Request, - db *sql.Conn, + conn *sql.Conn, ) (jr JSONResult, err error) { csi := input.(*models.CrossSectionInput) start := time.Now() - var rows *sql.Rows - if rows, err = db.QueryContext( - req.Context(), - crossSQL, - csi.Geometry.Coordinates.AsWKB(), - csi.Properties.Bottleneck, - csi.Properties.Date.Time, + tree, err := octree.Cache.Get( + csi.Properties.Bottleneck, csi.Properties.Date.Time, + conn, req.Context()) + + log.Printf("loading octree took: %s\n", time.Since(start)) + if err != nil { + return + } + + // The coordinate system of the octree is an UTM projection. + // The input coordinates are in WGS84. + // So we need to reproject them. + + start = time.Now() + + var rp *models.Reprojector + if rp, err = models.NewReprojector( + conn, req.Context(), + WGS84, tree.EPSG, ); err != nil { return } - defer rows.Close() - - var segments models.GeoJSONMultiLineCoordinatesZ + defer rp.Close() - for rows.Next() { - var segment models.GeoJSONLineCoordinatesZ - if err = rows.Scan(&segment); err != nil { - return - } - segments = append(segments, segment) - } + coords, err := reproject(rp, csi.Geometry.Coordinates, req.Context()) - if err = rows.Err(); err != nil { + log.Printf("transforming input coords took: %s\n", time.Since(start)) + if err != nil { return } - log.Printf("query took: %v\n", time.Since(start)) - start = time.Now() - joined := segments.Join() + var segments octree.MultiLineStringZ + + for i := 0; i < len(coords)-1; i++ { + c1 := &coords[i] + c2 := &coords[i+1] - log.Printf("joining took: %v\n", time.Since(start)) - log.Printf("segments before/after: %d %d\n", len(segments), len(joined)) + verticalLine := octree.NewVerticalLine(c1.Lat, c1.Lon, c2.Lat, c2.Lon) + + var line octree.MultiLineStringZ + tree.Vertical(c1.Lat, c1.Lon, c2.Lat, c2.Lon, func(t *octree.Triangle) { + if ls := verticalLine.Intersection(t); len(ls) > 0 { + line = append(line, ls) + } + }) - /* - if err2 := dumpProfile( - csi.Geometry.Coordinates[0], - csi.Geometry.Coordinates[len(csi.Geometry.Coordinates)-1], - joined, - ); err2 != nil { - log.Printf("error: %v\n", err2) + if len(line) > 0 { + log.Printf("line length: %d\n", len(line)) + // They are all on the segment (c1.Lat, c1.Lon) - (c2.Lat, c2.Lon). + // Sort them by project them on this line. + joined := line.JoinOnLine(c1.Lat, c1.Lon, c2.Lat, c2.Lon) + log.Printf("joined length: %d\n", len(joined)) + segments = append(segments, joined...) } - */ + + } + log.Printf("octree traversal took: %s\n", time.Since(start)) + + // The result have to be WGS84. So project the result back. + start = time.Now() + + rp.FromEPSG, rp.ToEPSG = tree.EPSG, WGS84 + + var joined models.GeoJSONMultiLineCoordinatesZ + joined, err = projectBack(rp, segments, req.Context()) + + log.Printf("projecting back took: %s\n", time.Since(start)) + if err != nil { + return + } jr = JSONResult{ Result: &models.CrossSectionOutput{ @@ -95,75 +150,3 @@ return } - -/* -func dumpProfile( - start models.GeoJSONCoordinate, - stop models.GeoJSONCoordinate, - segments models.GeoJSONMultiLineCoordinatesZ, -) error { - - w, err := os.Create("/tmp/cross.txt") - if err != nil { - return err - } - defer w.Close() - - begin := models.GeoJSONCoordinateZ{ - Lat: start.Lat, - Lon: start.Lon, - } - end := models.GeoJSONCoordinateZ{ - Lat: stop.Lat, - Lon: stop.Lon, - } - - last := begin - - var total float64 - - minz := 10000.0 - - for _, line := range segments { - for _, coord := range line { - if coord.Z < minz { - minz = coord.Z - } - total += last.Distance(coord) - last = coord - } - } - - log.Printf("length: %.3f minz: %f\n", total, minz) - - var pos float64 - - for i, line := range segments { - for j, coord := range line { - if i == 0 && j == 0 { - if !begin.Equals(coord) { - pos = begin.Distance(coord) - fmt.Fprintf(w, "%.3f %f\n", 0.0, 200.0) - fmt.Fprintf(w, "%.3f %f\n", pos-0.0001, 200.0) - fmt.Fprintf(w, "%.3f %f\n", pos, coord.Z-minz) - continue - } - } else if j == 0 { - fmt.Fprintf(w, "%.3f %f\n", pos+0.0001, 200.0) - fmt.Fprintf(w, "%.3f %f\n", pos+last.Distance(coord)-0.0001, 200.0) - continue - } - pos += last.Distance(coord) - fmt.Fprintf(w, "%.3f %f\n", pos, coord.Z-minz) - last = coord - } - } - - if !last.Equals(end) { - fmt.Fprintf(w, "%.3f %f\n", pos+0.0001, 200.0) - fmt.Fprintf(w, "%.3f %f\n", pos+last.Distance(end), 200.0) - } - - return nil -} -*/ diff -r 0689f4b7c99f -r b6a1779ffb42 pkg/controllers/octreecross.go --- a/pkg/controllers/octreecross.go Thu Sep 27 12:49:45 2018 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,152 +0,0 @@ -package controllers - -import ( - "context" - "database/sql" - "log" - "net/http" - "time" - - "gemma.intevation.de/gemma/pkg/models" - "gemma.intevation.de/gemma/pkg/octree" -) - -const WGS84 = 4326 - -func reproject( - rp *models.Reprojector, - src models.GeoJSONLineCoordinates, - ctx context.Context, -) (models.GeoJSONLineCoordinates, error) { - - dst := make(models.GeoJSONLineCoordinates, len(src)) - for i, s := range src { - var err error - if dst[i].Lat, dst[i].Lon, err = rp.Reproject( - s.Lat, s.Lon, ctx, - ); err != nil { - return nil, err - } - } - return dst, nil -} - -func projectBack( - rp *models.Reprojector, - lines octree.MultiLineStringZ, - ctx context.Context, -) (models.GeoJSONMultiLineCoordinatesZ, error) { - - out := make(models.GeoJSONMultiLineCoordinatesZ, len(lines)) - - for i, segment := range lines { - - coords := make(models.GeoJSONLineCoordinatesZ, len(segment)) - out[i] = coords - for j, v := range segment { - lat, lon, err := rp.Reproject(v.X, v.Y, ctx) - if err != nil { - return nil, err - } - coords[j] = models.GeoJSONCoordinateZ{Lat: lat, Lon: lon, Z: v.Z} - } - } - return out, nil -} - -func octreeCrossSection( - input interface{}, - req *http.Request, - conn *sql.Conn, -) (jr JSONResult, err error) { - - csi := input.(*models.CrossSectionInput) - - start := time.Now() - - tree, err := octree.Cache.Get( - csi.Properties.Bottleneck, csi.Properties.Date.Time, - conn, req.Context()) - - log.Printf("loading octree took: %s\n", time.Since(start)) - if err != nil { - return - } - - // The coordinate system of the octree is an UTM projection. - // The input coordinates are in WGS84. - // So we need to reproject them. - - start = time.Now() - - var rp *models.Reprojector - if rp, err = models.NewReprojector( - conn, req.Context(), - WGS84, tree.EPSG, - ); err != nil { - return - } - defer rp.Close() - - coords, err := reproject(rp, csi.Geometry.Coordinates, req.Context()) - - log.Printf("transforming input coords took: %s\n", time.Since(start)) - if err != nil { - return - } - - start = time.Now() - - var segments octree.MultiLineStringZ - - for i := 0; i < len(coords)-1; i++ { - c1 := &coords[i] - c2 := &coords[i+1] - - verticalLine := octree.NewVerticalLine(c1.Lat, c1.Lon, c2.Lat, c2.Lon) - - var line octree.MultiLineStringZ - tree.Vertical(c1.Lat, c1.Lon, c2.Lat, c2.Lon, func(t *octree.Triangle) { - if ls := verticalLine.Intersection(t); len(ls) > 0 { - line = append(line, ls) - } - }) - - if len(line) > 0 { - log.Printf("line length: %d\n", len(line)) - // They are all on the segment (c1.Lat, c1.Lon) - (c2.Lat, c2.Lon). - // Sort them by project them on this line. - joined := line.JoinOnLine(c1.Lat, c1.Lon, c2.Lat, c2.Lon) - log.Printf("joined length: %d\n", len(joined)) - segments = append(segments, joined...) - } - - } - log.Printf("octree traversal took: %s\n", time.Since(start)) - - // The result have to be WGS84. So project the result back. - start = time.Now() - - rp.FromEPSG, rp.ToEPSG = tree.EPSG, WGS84 - - var joined models.GeoJSONMultiLineCoordinatesZ - joined, err = projectBack(rp, segments, req.Context()) - - log.Printf("projecting back took: %s\n", time.Since(start)) - if err != nil { - return - } - - jr = JSONResult{ - Result: &models.CrossSectionOutput{ - Type: "Feature", - Geometry: models.CrossSectionOutputGeometry{ - Type: "MultiLineString", - Coordinates: joined, - }, - Properties: map[string]interface{}{}, - }, - } - - return -} diff -r 0689f4b7c99f -r b6a1779ffb42 pkg/controllers/routes.go --- a/pkg/controllers/routes.go Thu Sep 27 12:49:45 2018 +0200 +++ b/pkg/controllers/routes.go Thu Sep 27 13:18:50 2018 +0200 @@ -105,11 +105,6 @@ // Cross sections - api.Handle("/octcross", any(&JSONHandler{ - Input: func() interface{} { return new(models.CrossSectionInput) }, - Handle: octreeCrossSection, - })).Methods(http.MethodPost) - api.Handle("/cross", any(&JSONHandler{ Input: func() interface{} { return new(models.CrossSectionInput) }, Handle: crossSection, diff -r 0689f4b7c99f -r b6a1779ffb42 pkg/models/cross.go --- a/pkg/models/cross.go Thu Sep 27 12:49:45 2018 +0200 +++ b/pkg/models/cross.go Thu Sep 27 13:18:50 2018 +0200 @@ -7,7 +7,6 @@ "errors" "fmt" "math" - "sort" "time" ) @@ -257,121 +256,3 @@ x := dLng * math.Cos(deg2rad((cz.Lat+other.Lat)/2.0)) return math.Sqrt(dLat*dLat+x*x) * EarthRadius } - -type point struct { - x float64 - y float64 -} - -func (cz GeoJSONCoordinateZ) quant() point { - const xyScale = 1e6 - return point{ - math.Round(cz.Lon*xyScale) / xyScale, - math.Round(cz.Lat*xyScale) / xyScale, - } -} - -func (lcz GeoJSONLineCoordinatesZ) join(other GeoJSONLineCoordinatesZ) GeoJSONLineCoordinatesZ { - l1, l2 := len(lcz), len(other) - n := make(GeoJSONLineCoordinatesZ, l1+l2-1) - copy(n, lcz[:l1-1]) - n[l1-1] = lcz[l1-1].mid(other[0]) - copy(n[l1:], other[1:]) - return n -} - -func (ml GeoJSONMultiLineCoordinatesZ) Join() GeoJSONMultiLineCoordinatesZ { - - if len(ml) < 2 { - return ml - } - - type value struct { - coords GeoJSONLineCoordinatesZ - order int - } - - heads := make(map[point]*value) - - for i, coords := range ml { - key := coords[len(coords)-1].quant() - if v := heads[key]; v != nil { - delete(heads, key) - v.coords = coords.join(v.coords) - heads[v.coords[0].quant()] = v - } else { - heads[coords[0].quant()] = &value{coords, i} - } - } - -again: - for k, v1 := range heads { - key := v1.coords[len(v1.coords)-1].quant() - if k == key { // cycle detected - continue - } - if v2 := heads[key]; v2 != nil { - delete(heads, key) - v1.coords = v1.coords.join(v2.coords) - goto again - } - } - - // To make it deterministic sort in input order. - tmp := make([]*value, len(heads)) - var i int - for _, v := range heads { - tmp[i] = v - i++ - } - sort.Slice(tmp, func(i, j int) bool { return tmp[i].order < tmp[j].order }) - - out := make(GeoJSONMultiLineCoordinatesZ, len(tmp)) - for i, v := range tmp { - out[i] = v.coords - } - - return out -} - -/* -func (ml GeoJSONMultiLineCoordinatesZ) Join() GeoJSONMultiLineCoordinatesZ { - - if len(ml) < 2 { - return ml - } - - lists := make(GeoJSONMultiLineCoordinatesZ, len(ml)) - copy(lists, ml) - -next: - for j := 0; j < len(lists); { - front := lists[j] - tail := front[len(front)-1] - - for i := 0; i < len(lists); i++ { - if i == j { - continue - } - curr := lists[i] - head := curr[0] - if tail.Equals(head) { - n := make(GeoJSONLineCoordinatesZ, len(front)+len(curr)-1) - copy(n, front[:len(front)-1]) - n[len(front)-1] = head.mid(tail) - copy(n[len(front):], curr[1:]) - lists[j] = n - if i < len(lists)-1 { - copy(lists[i:], lists[i+1:]) - } - lists[len(lists)-1] = nil - lists = lists[:len(lists)-1] - continue next - } - } - j++ - } - - return lists -} -*/