Merge pull request #47 from Dewberry/feature/mesh-extraction

Feature/mesh extraction

Dockerfile

@@ -34,7 +34,6 @@ RUN rm HEC-RAS_62_Example_Projects.zip
 
 ENTRYPOINT /app/main
 
-
 FROM ghcr.io/osgeo/gdal:alpine-normal-3.8.3 as prod
 COPY --from=build /app/main /main
 ENTRYPOINT /main

LICENSE

@@ -2,6 +2,10 @@ MIT License
 
 Copyright (c) 2020 USACE
 
+Copyright (C) 2013 Przemyslaw Szczepaniak https://github.com/pzsz/voronoi
+
+Copyright (c) 2010-2013 Raymond Hill https://github.com/gorhill/Javascript-Voronoi
+
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal
 in the Software without restriction, including without limitation the rights

go.mod

@@ -10,6 +10,7 @@ require (
 	github.com/jackc/pgx v3.6.2+incompatible
 	github.com/jmoiron/sqlx v1.3.4
 	github.com/labstack/echo/v4 v4.3.0
+	github.com/pzsz/voronoi v0.0.0-20130609164533-4314be88c79f
 	github.com/swaggo/swag v1.7.0
 )
 

go.sum

@@ -90,6 +90,8 @@ github.com/pkg/errors v0.9.1 h1:FEBLx1zS214owpjy7qsBeixbURkuhQAwrK5UwLGTwt4=
 github.com/pkg/errors v0.9.1/go.mod h1:bwawxfHBFNV+L2hUp1rHADufV3IMtnDRdf1r5NINEl0=
 github.com/pmezard/go-difflib v1.0.0 h1:4DBwDE0NGyQoBHbLQYPwSUPoCMWR5BEzIk/f1lZbAQM=
 github.com/pmezard/go-difflib v1.0.0/go.mod h1:iKH77koFhYxTK1pcRnkKkqfTogsbg7gZNVY4sRDYZ/4=
+github.com/pzsz/voronoi v0.0.0-20130609164533-4314be88c79f h1:YeUeWAPTrycaikVXiVoAu1dxIR4+tSsn9IbTkNfRsss=
+github.com/pzsz/voronoi v0.0.0-20130609164533-4314be88c79f/go.mod h1:T6EcHUIQ7r0Xr1jjOJ7hmpUaqVyqb9WimeNzS9nz4Zo=
 github.com/russross/blackfriday/v2 v2.0.1/go.mod h1:+Rmxgy9KzJVeS9/2gXHxylqXiyQDYRxCVz55jmeOWTM=
 github.com/shopspring/decimal v1.2.0 h1:abSATXmQEYyShuxI4/vyW3tV1MrKAJzCZ/0zLUXYbsQ=
 github.com/shopspring/decimal v1.2.0/go.mod h1:DKyhrW/HYNuLGql+MJL6WCR6knT2jwCFRcu2hWCYk4o=

tools/geospatial.go

@@ -12,6 +12,7 @@ import (
 	"github.com/USACE/filestore"
 	"github.com/dewberry/gdal"
 	"github.com/go-errors/errors" // warning: replaces standard errors
+	"github.com/pzsz/voronoi"
 )
 
 // GeoData ...
@@ -27,6 +28,7 @@ type Features struct {
 	Banks               []VectorFeature
 	StorageAreas        []VectorFeature
 	TwoDAreas           []VectorFeature
+	Mesh                []VectorFeature
 	HydraulicStructures []VectorFeature
 	Connections         []VectorFeature
 	BCLines             []VectorFeature
@@ -206,7 +208,7 @@ func getRiverCenterline(sc *bufio.Scanner, transform gdal.CoordinateTransform) (
 	}
 
 	xyLineString.Transform(transform)
-	// This is a temporary fix since the x and y values need to be flipped:
+	// The x and y values need to be flipped:
 	yxLineString := flipXYLineString(xyLineString)
 
 	multiLineString := yxLineString.ForceToMultiLineString()
@@ -291,7 +293,7 @@ func getXS(sc *bufio.Scanner, transform gdal.CoordinateTransform, riverReachName
 	}
 
 	xyzLineString.Transform(transform)
-	// This is a temporary fix since the x and y values need to be flipped
+	// The x and y values need to be flipped
 	yxzLineString := flipXYLineString25D(xyzLineString)
 
 	multiLineString := yxzLineString.ForceToMultiLineString()
@@ -319,7 +321,7 @@ func getBanks(line string, transform gdal.CoordinateTransform, xsFeature VectorF
 		xyPoint := gdal.Create(gdal.GT_Point)
 		xyPoint.AddPoint2D(bankXY[0], bankXY[1])
 		xyPoint.Transform(transform)
-		// This is a temporary fix since the x and y values need to be flipped
+		// The x and y values need to be flipped
 		yxPoint := flipXYPoint(xyPoint)
 		multiPoint := yxPoint.ForceToMultiPoint()
 		wkb, err := multiPoint.ToWKB()
@@ -351,7 +353,7 @@ func getArea(sc *bufio.Scanner, transform gdal.CoordinateTransform) (VectorFeatu
 	}
 
 	xyLinearRing.Transform(transform)
-	// This is a temporary fix since the x and y values need to be flipped:
+	// The x and y values need to be flipped
 	yxLinearRing := flipXYLinearRing(xyLinearRing)
 
 	yxPolygon := gdal.Create(gdal.GT_Polygon)
@@ -365,6 +367,75 @@ func getArea(sc *bufio.Scanner, transform gdal.CoordinateTransform) (VectorFeatu
 	return feature, is2D, nil
 }
 
+func validVoronoiBound(point [2]float64, bbox voronoi.BBox) bool {
+	return point[0] > bbox.Xl && point[0] < bbox.Xr && point[1] > bbox.Yt && point[1] < bbox.Yb
+}
+
+func getMeshArea(sc *bufio.Scanner, transform gdal.CoordinateTransform, allowedDist float64) ([]VectorFeature, error) {
+	features := []VectorFeature{
+		VectorFeature{FeatureName: "mesh_points"},
+		VectorFeature{FeatureName: "mesh_voronoi"},
+		VectorFeature{FeatureName: "mesh_concave"},
+	}
+
+	xyPairs, err := getDataPairsfromTextBlock("Storage Area 2D Points=", sc, 64, 16)
+	if err != nil {
+		return features, errors.Wrap(err, 0)
+	}
+
+	multipoint := gdal.Create(gdal.GT_MultiPoint) // for multipoint
+	vertices := voronoi.Vertices{}                // for voronoi
+
+	for _, point := range xyPairs {
+		// Multipoint creation
+		xyPoint := gdal.Create(gdal.GT_Point)
+		xyPoint.AddPoint2D(point[0], point[1])
+		xyPoint.Transform(transform)
+		// The x and y values need to be flipped
+		yxPoint := flipXYPoint(xyPoint)
+
+		err = multipoint.AddGeometry(yxPoint)
+		if err != nil {
+			return features, errors.Wrap(err, 0)
+		}
+
+		// Voronoi Vertex population
+		vertex := voronoi.Vertex{X: yxPoint.X(0), Y: yxPoint.Y(0)}
+		vertices = append(vertices, vertex)
+	}
+
+	concaveHull := multipoint.ConcaveHull(.02, false)
+	bounds := concaveHull.Envelope()
+
+	bbox := voronoi.NewBBox(bounds.MinX(), bounds.MaxX(), bounds.MinY(), bounds.MaxY())
+	diagram := voronoi.ComputeDiagram(vertices, bbox, false)
+
+	voronoiMultiLineString := gdal.Create(gdal.GT_MultiLineString)
+	for _, edge := range diagram.Edges {
+		lineString := gdal.Create(gdal.GT_LineString)
+		vStart := edge.Va.Vertex
+		vEnd := edge.Vb.Vertex
+
+		lineString.AddPoint2D(vStart.X, vStart.Y)
+		lineString.AddPoint2D(vEnd.X, vEnd.Y)
+
+		if concaveHull.Contains(lineString) {
+			voronoiMultiLineString.AddGeometry(lineString)
+		}
+		lineString.Destroy()
+	}
+
+	for ind, geom := range []gdal.Geometry{multipoint, voronoiMultiLineString, concaveHull} {
+		wkb, err := geom.ToWKB()
+		if err != nil {
+			return features, errors.Wrap(err, 0)
+		}
+		features[ind].Geometry = wkb
+	}
+
+	return features, nil
+}
+
 func getAreaType(sc *bufio.Scanner) (string, error) {
 	is2D := ""
 
@@ -403,7 +474,7 @@ func getBreakLine(sc *bufio.Scanner, transform gdal.CoordinateTransform) (Vector
 	}
 
 	xyLineString.Transform(transform)
-	// This is a temporary fix since the x and y values need to be flipped:
+	// The x and y values need to be flipped
 	yxLineString := flipXYLineString(xyLineString)
 
 	multiLineString := yxLineString.ForceToMultiLineString()
@@ -447,7 +518,7 @@ func getBCLine(sc *bufio.Scanner, transform gdal.CoordinateTransform) (VectorFea
 	}
 
 	xyLineString.Transform(transform)
-	// This is a temporary fix since the x and y values need to be flipped:
+	// The x and y values need to be flipped
 	yxLineString := flipXYLineString(xyLineString)
 
 	multiLineString := yxLineString.ForceToMultiLineString()
@@ -496,7 +567,7 @@ func getConnectionLine(sc *bufio.Scanner, transform gdal.CoordinateTransform) (V
 	}
 
 	xyLineString.Transform(transform)
-	// This is a temporary fix since the x and y values need to be flipped:
+	// The x and y values need to be flipped
 	yxLineString := flipXYLineString(xyLineString)
 
 	multiLineString := yxLineString.ForceToMultiLineString()
@@ -577,6 +648,14 @@ func GetGeospatialData(gd *GeoData, fs filestore.FileStore, geomFilePath string,
 			} else if aType == "-1" {
 				f.TwoDAreas = append(f.TwoDAreas, storageAreaFeature)
 			}
+		case strings.HasPrefix(line, "Storage Area 2D Points="):
+			epsilon := 1e-1
+			meshFeatures, err := getMeshArea(sc, transform, epsilon)
+			if err != nil {
+				return errors.Wrap(err, 0)
+			}
+			f.Mesh = append(f.Mesh, meshFeatures...)
+
 		case strings.HasPrefix(line, "Type RM Length L Ch R = 1"):
 			xsFeature, bankLayer, err := getXSBanks(sc, transform, riverReachName)
 			if err != nil {