diff --git a/config/config.go b/config/config.go index b61f6a1..9199692 100644 --- a/config/config.go +++ b/config/config.go @@ -8,9 +8,10 @@ import ( ) type APIConfig struct { - Host string - Port int - FileStore *filestore.FileStore + Host string + Port int + FileStore *filestore.FileStore + DestinationCRS int } // Address tells the application where to run the api out of @@ -24,6 +25,7 @@ func Init() *APIConfig { config.Host = "" // 0.0.0.0 config.Port = 5600 config.FileStore = FileStoreInit(false) + config.DestinationCRS = 4326 return config } diff --git a/handlers/geospatialdata.go b/handlers/geospatialdata.go index 8bb581d..3aeada3 100644 --- a/handlers/geospatialdata.go +++ b/handlers/geospatialdata.go @@ -3,9 +3,9 @@ package handlers import ( "net/http" + "github.com/USACE/mcat-ras/config" ras "github.com/USACE/mcat-ras/tools" - "github.com/USACE/filestore" "github.com/labstack/echo/v4" ) @@ -19,17 +19,17 @@ import ( // @Success 200 {object} interface{} // @Failure 500 {object} SimpleResponse // @Router /geospatialdata [get] -func GeospatialData(fs *filestore.FileStore) echo.HandlerFunc { +func GeospatialData(ac *config.APIConfig) echo.HandlerFunc { return func(c echo.Context) error { definitionFile := c.QueryParam("definition_file") - rm, err := ras.NewRasModel(definitionFile, *fs) + rm, err := ras.NewRasModel(definitionFile, *ac.FileStore) if err != nil { return c.JSON(http.StatusInternalServerError, SimpleResponse{http.StatusInternalServerError, err.Error()}) } - data, err := rm.GeospatialData() + data, err := rm.GeospatialData(ac.DestinationCRS) if err != nil { return c.JSON(http.StatusInternalServerError, SimpleResponse{http.StatusInternalServerError, err.Error()}) } diff --git a/handlers/index.go b/handlers/index.go index e994be6..79e3f2e 100644 --- a/handlers/index.go +++ b/handlers/index.go @@ -28,11 +28,11 @@ func Index(fs *filestore.FileStore) echo.HandlerFunc { if err != nil { return c.JSON(http.StatusInternalServerError, SimpleResponse{http.StatusInternalServerError, err.Error()}) } - mod, err := rm.Index() - if err != nil { - return c.JSON(http.StatusInternalServerError, SimpleResponse{http.StatusInternalServerError, err.Error()}) - } + // mod, err := rm.Index() + // if err != nil { + // return c.JSON(http.StatusInternalServerError, SimpleResponse{http.StatusInternalServerError, err.Error()}) + // } - return c.JSON(http.StatusOK, mod) + return c.JSON(http.StatusOK, rm) } } diff --git a/main.go b/main.go index 743bc50..d413225 100644 --- a/main.go +++ b/main.go @@ -39,7 +39,7 @@ func main() { e.GET("/isgeospatial", handlers.IsGeospatial(appConfig.FileStore)) e.GET("/modeltype", handlers.ModelType(appConfig.FileStore)) e.GET("/modelversion", handlers.ModelVersion(appConfig.FileStore)) - e.GET("/geospatialdata", handlers.GeospatialData(appConfig.FileStore)) + e.GET("/geospatialdata", handlers.GeospatialData(appConfig)) e.Logger.Fatal(e.Start(appConfig.Address())) } diff --git a/tools/geom.go b/tools/geom.go index bc3b20e..9865173 100644 --- a/tools/geom.go +++ b/tools/geom.go @@ -2,26 +2,20 @@ package tools import ( "bufio" - "errors" "fmt" - "math" "path/filepath" - "regexp" - "strconv" "strings" "sync" - - "github.com/USACE/filestore" - "github.com/dewberry/gdal" ) // GeomFileContents keywords and data container for ras flow file search type GeomFileContents struct { Path string - FileExt string //`json:"File Extension"` - GeomTitle string //`json:"Geom Title"` - ProgramVersion string //`json:"Program Version"` - Description string //`json:"Description"` + FileExt string `json:"File Extension"` + GeomTitle string `json:"Geom Title"` + ProgramVersion string `json:"Program Version"` + Description string `json:"Description"` + Structures []hydraulicStructures `json:"Hydraulic Structures"` } // getGeomData Reads a geometry file. returns none to allow concurrency @@ -39,414 +33,44 @@ func getGeomData(rm *RasModel, fn string, wg *sync.WaitGroup, errChan chan error defer f.Close() sc := bufio.NewScanner(f) - var line string - header := true - for sc.Scan() { - - line = sc.Text() - - match, err := regexp.MatchString("=", line) - - if err != nil { - errChan <- err - return - } - - beginDescription, err := regexp.MatchString("BEGIN GEOM DESCRIPTION", line) - - if err != nil { - errChan <- err - return - } - - if match { - data := strings.Split(line, "=") - - switch data[0] { - - case "Geom Title": - meta.GeomTitle = data[1] - - case "Program Version": - meta.ProgramVersion = data[1] - - case "River Reach", "Storage Area": - header = false - - } - - } else if header && beginDescription { - - for sc.Scan() { - line = sc.Text() - endDescription, _ := regexp.MatchString("END GEOM DESCRIPTION", line) - - if endDescription { - break - - } else { - if line != "" { - meta.Description += line + "\n" - } - } - - } - - } - } - - rm.Metadata.GeomFiles = append(rm.Metadata.GeomFiles, meta) - return -} - -// sourceCRS -var sourceCRS string = `PROJCS["NAD_1983_StatePlane_Texas_Central_FIPS_4203_Feet",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",2296583.333333333],PARAMETER["False_Northing",9842500.0],PARAMETER["Central_Meridian",-100.3333333333333],PARAMETER["Standard_Parallel_1",30.11666666666667],PARAMETER["Standard_Parallel_2",31.88333333333333],PARAMETER["Latitude_Of_Origin",29.66666666666667],UNIT["Foot_US",0.3048006096012192]]` - -// DestinationCRS ... -var DestinationCRS int = 4326 - -// GeoData ... -type GeoData struct { - Features map[string]Features - Georeference int -} - -// Features ... -type Features struct { - Rivers []vectorLayer - XS []vectorLayer - Banks []vectorLayer - StorageAreas []vectorLayer - TwoDAreas []vectorLayer - HydraulicStructures []vectorLayer -} - -type vectorLayer struct { - FeatureName string `json:"feature_name"` - Fields map[string]interface{} `json:"fields"` - Geometry []uint8 `json:"geometry"` -} -type xyzPoint struct { - x float64 - y float64 - z float64 -} + var description string -func rightofEquals(line string) string { - return strings.TrimSpace(strings.Split(line, "=")[1]) -} - -func dataPairsfromTextBlock(sc *bufio.Scanner, nPairs int, colWidth int, valueWidth int) ([][2]float64, error) { - var stride int = valueWidth * 2 - pairs := [][2]float64{} -out: - for sc.Scan() { - line := sc.Text() - for s := 0; s < colWidth; { - if len(line) > s { - val1, err := strconv.ParseFloat(strings.TrimSpace(line[s:s+valueWidth]), 64) - if err != nil { - return pairs, err - } - val2, err := strconv.ParseFloat(strings.TrimSpace(line[s+valueWidth:s+stride]), 64) - if err != nil { - return pairs, err - } - pairs = append(pairs, [2]float64{val1, val2}) - if int(len(pairs)) == nPairs { - break out - } - } else { - break - } - s += stride - } - } - return pairs, nil -} - -func getDataPairsfromTextBlock(nDataPairsLine string, sc *bufio.Scanner, colWidth int, valueWidth int) ([][2]float64, error) { - pairs := [][2]float64{} + header := true + idx := 0 for sc.Scan() { + idx++ line := sc.Text() - if strings.HasPrefix(line, nDataPairsLine) { - nPairs, err := strconv.Atoi(rightofEquals(line)) - if err != nil { - return pairs, err - } - pairs, err = dataPairsfromTextBlock(sc, nPairs, colWidth, valueWidth) - if err != nil { - return pairs, err - } - break - } - } - return pairs, nil -} - -// distance returns the distance along a straight line in euclidean space -func distance(p0, p1 [2]float64) float64 { - result := math.Sqrt(math.Pow((p1[0]-p0[0]), 2) + math.Pow((p1[1]-p0[1]), 2)) - return result -} - -// pointAtDistance returns a new point along a straight line in euclidean space -// at a specified distance -func pointAtDistance(p0, p1 [2]float64, delta float64) [2]float64 { - distanceRatio := delta / distance(p0, p1) - newX := (1-distanceRatio)*p0[0] + distanceRatio*p1[0] - newY := (1-distanceRatio)*p0[1] + distanceRatio*p1[1] - return [2]float64{newX, newY} -} - -// interpZ creates a new point a given distance along a line composed -// of many segments. -func interpXY(xyPairs [][2]float64, d float64) [2]float64 { - // newPoint is an x, y pair - var newPoint [2]float64 - lineSegments := len(xyPairs) - 1 - lineLength := 0.0 - -findLineSegment: - for i := 0; i < lineSegments; i++ { - p0, p1 := xyPairs[i], xyPairs[i+1] - lineLength += distance(p0, p1) - switch { - case lineLength > d: - delta := distance(p0, p1) - (lineLength - d) - newPoint = pointAtDistance(p0, p1, delta) - break findLineSegment + case strings.HasPrefix(line, "Geom Title="): + meta.GeomTitle = rightofEquals(line) - default: - continue - } - } - if d >= lineLength { - if d-lineLength <= 0.1 { - return xyPairs[len(xyPairs)-1] - } - fmt.Printf("The interpolated point has a station of %v while the xy line is %v long", d, lineLength) - } - return newPoint -} + case strings.HasPrefix(line, "Program Version="): + meta.ProgramVersion = rightofEquals(line) -// attributeZ using station from cross-section line and gis coordinates -func attributeZ(xyPairs [][2]float64, mzPairs [][2]float64) []xyzPoint { - points := []xyzPoint{} - startingStation := mzPairs[0][0] - for _, mzPair := range mzPairs { - newPoint := interpXY(xyPairs, mzPair[0]-startingStation) - if newPoint[0] != 0 && newPoint[1] != 0 { - points = append(points, xyzPoint{newPoint[0], newPoint[1], mzPair[1]}) - } else { - fmt.Printf("Interpolated point has an xy value of (%v, %v). ", newPoint[0], newPoint[1]) - } - } - return points -} - -func getTransform(sourceProjection string, destinationEPSG int) (gdal.CoordinateTransform, error) { - transform := gdal.CoordinateTransform{} - sourceSpRef := gdal.CreateSpatialReference(sourceProjection) - sourceSpRef.MorphFromESRI() - if err := sourceSpRef.Validate(); err != nil { - return transform, errors.New("Unable to extract geospatial data. Invalid source Projection") - } - - destinationSpRef := gdal.CreateSpatialReference("") - if err := destinationSpRef.FromEPSG(destinationEPSG); err != nil { - return transform, err - } - transform = gdal.CreateCoordinateTransform(sourceSpRef, destinationSpRef) - return transform, nil -} - -func flipXYLineString(xyLineString gdal.Geometry) gdal.Geometry { - yxLineString := gdal.Create(gdal.GT_LineString) - nPoints := xyLineString.PointCount() - for i := 0; i < nPoints; i++ { - x, y, _ := xyLineString.Point(i) - yxLineString.AddPoint2D(y, x) - } - xyLineString.Destroy() - return yxLineString -} - -func flipXYLineString25D(xyzLineString gdal.Geometry) gdal.Geometry { - yxzLineString := gdal.Create(gdal.GT_LineString25D) - nPoints := xyzLineString.PointCount() - for i := 0; i < nPoints; i++ { - x, y, z := xyzLineString.Point(i) - yxzLineString.AddPoint(y, x, z) - } - xyzLineString.Destroy() - return yxzLineString -} - -func flipXYPoint(xyPoint gdal.Geometry) gdal.Geometry { - yxPoint := gdal.Create(gdal.GT_Point) - nPoints := xyPoint.PointCount() - for i := 0; i < nPoints; i++ { - x, y, _ := xyPoint.Point(i) - yxPoint.AddPoint2D(y, x) - } - xyPoint.Destroy() - return yxPoint -} - -func getRiverCenterline(sc *bufio.Scanner, transform gdal.CoordinateTransform) (vectorLayer, error) { - riverReach := strings.Split(rightofEquals(sc.Text()), ",") - layer := vectorLayer{FeatureName: fmt.Sprintf("%s, %s", strings.TrimSpace(riverReach[0]), strings.TrimSpace(riverReach[1]))} - - xyPairs, err := getDataPairsfromTextBlock("Reach XY=", sc, 64, 16) - if err != nil { - return layer, err - } - - xyLineString := gdal.Create(gdal.GT_LineString) - for _, pair := range xyPairs { - xyLineString.AddPoint2D(pair[0], pair[1]) - } - - xyLineString.Transform(transform) - // This is a temporary fix since the x and y values need to be flipped: - yxLineString := flipXYLineString(xyLineString) - - multiLineString := yxLineString.ForceToMultiLineString() - wkb, err := multiLineString.ToWKB() - if err != nil { - return layer, err - } - layer.Geometry = wkb - return layer, err -} - -func getXSBanks(sc *bufio.Scanner, transform gdal.CoordinateTransform, riverReachName string) (vectorLayer, []vectorLayer, error) { - bankLayers := []vectorLayer{} - - xsLayer, xyPairs, startingStation, err := getXS(sc, transform, riverReachName) - if err != nil { - return xsLayer, bankLayers, err - } - for sc.Scan() { - line := sc.Text() - if strings.HasPrefix(line, "Bank Sta=") { - bankLayers, err = getBanks(line, transform, xsLayer, xyPairs, startingStation) - if err != nil { - return xsLayer, bankLayers, err + case strings.HasPrefix(line, "BEGIN GEOM DESCRIPTION:"): + if header { + description, idx, err = getDescription(sc, idx, "END GEOM DESCRIPTION:") + if err != nil { + errChan <- err + return + } + meta.Description += description } - break - } - } - return xsLayer, bankLayers, err -} -func getXS(sc *bufio.Scanner, transform gdal.CoordinateTransform, riverReachName string) (vectorLayer, [][2]float64, float64, error) { - compData := strings.Split(rightofEquals(sc.Text()), ",") - layer := vectorLayer{FeatureName: strings.TrimSpace(compData[1]), Fields: map[string]interface{}{}} - layer.Fields["RiverReachName"] = riverReachName - - xyPairs, err := getDataPairsfromTextBlock("XS GIS Cut Line", sc, 64, 16) - if err != nil { - return layer, xyPairs, 0.0, err - } - - mzPairs, err := getDataPairsfromTextBlock("#Sta/Elev", sc, 80, 8) - if err != nil { - return layer, xyPairs, mzPairs[0][0], err - } - - xyzPoints := attributeZ(xyPairs, mzPairs) - - xyzLineString := gdal.Create(gdal.GT_LineString25D) - for _, point := range xyzPoints { - xyzLineString.AddPoint(point.x, point.y, point.z) - } - - xyzLineString.Transform(transform) - // This is a temporary fix since the x and y values need to be flipped - yxzLineString := flipXYLineString25D(xyzLineString) - - multiLineString := yxzLineString.ForceToMultiLineString() - wkb, err := multiLineString.ToWKB() - if err != nil { - return layer, xyPairs, mzPairs[0][0], err - } - layer.Geometry = wkb - return layer, xyPairs, mzPairs[0][0], err -} - -func getBanks(line string, transform gdal.CoordinateTransform, xsLayer vectorLayer, xyPairs [][2]float64, startingStation float64) ([]vectorLayer, error) { - layers := []vectorLayer{} - - bankStations := strings.Split(rightofEquals(line), ",") - for _, s := range bankStations { - layer := vectorLayer{FeatureName: strings.TrimSpace(s), Fields: map[string]interface{}{}} - layer.Fields["RiverReachName"] = xsLayer.Fields["RiverReachName"] - layer.Fields["xsName"] = xsLayer.FeatureName - bankStation, err := strconv.ParseFloat(s, 64) - if err != nil { - return layers, err - } - bankXY := interpXY(xyPairs, bankStation-startingStation) - 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 - yxPoint := flipXYPoint(xyPoint) - multiPoint := yxPoint.ForceToMultiPoint() - wkb, err := multiPoint.ToWKB() - if err != nil { - return layers, err - } - layer.Geometry = wkb - layers = append(layers, layer) - } - return layers, nil -} - -// GetGeospatialData ... -func GetGeospatialData(gd *GeoData, fs filestore.FileStore, geomFilePath string) error { - geomFileName := filepath.Base(geomFilePath) - f := Features{} - riverReachName := "" - - file, err := fs.GetObject(geomFilePath) - if err != nil { - return err - } - defer file.Close() - - sc := bufio.NewScanner(file) - - transform, err := getTransform(sourceCRS, DestinationCRS) - if err != nil { - return err - } - - for sc.Scan() { - line := sc.Text() - switch { case strings.HasPrefix(line, "River Reach="): - riverLayer, err := getRiverCenterline(sc, transform) + structures, err := getHydraulicStructureData(rm, fn, idx) if err != nil { - return err + fmt.Println(err) + continue } - f.Rivers = append(f.Rivers, riverLayer) - riverReachName = riverLayer.FeatureName + meta.Structures = append(meta.Structures, structures) + header = false - case strings.HasPrefix(line, "Type RM Length L Ch R = 1"): - xsLayer, bankLayers, err := getXSBanks(sc, transform, riverReachName) - if err != nil { - return err - } - f.XS = append(f.XS, xsLayer) - f.Banks = append(f.Banks, bankLayers...) + case strings.HasPrefix(line, "Storage Area="): + header = false } } - - gd.Features[geomFileName] = f - return nil + rm.Metadata.GeomFiles = append(rm.Metadata.GeomFiles, meta) + return } diff --git a/tools/geospatial.go b/tools/geospatial.go new file mode 100644 index 0000000..143410e --- /dev/null +++ b/tools/geospatial.go @@ -0,0 +1,416 @@ +package tools + +import ( + "bufio" + "errors" + "fmt" + "math" + "path/filepath" + "strconv" + "strings" + + "github.com/USACE/filestore" + "github.com/dewberry/gdal" +) + +// GeoData ... +type GeoData struct { + Features map[string]Features + Georeference int +} + +// Features ... +type Features struct { + Rivers []VectorLayer + XS []VectorLayer + Banks []VectorLayer + StorageAreas []VectorLayer + TwoDAreas []VectorLayer + HydraulicStructures []VectorLayer +} + +// VectorLayer ... +type VectorLayer struct { + FeatureName string `json:"feature_name"` + Fields map[string]interface{} `json:"fields"` + Geometry []uint8 `json:"geometry"` +} + +type xyzPoint struct { + x float64 + y float64 + z float64 +} + +var unitConsistencyMap map[string]string = map[string]string{ + "English Units": "US survey foot", + "SI Units": "metre"} + +// checkUnitConsistency checks that the unit system used by the model and its coordinate reference system are the same +func checkUnitConsistency(modelUnits string, sourceCRS string) error { + sourceSpRef := gdal.CreateSpatialReference(sourceCRS) + + if crsUnits, ok := sourceSpRef.AttrValue("UNIT", 0); ok { + if unitConsistencyMap[modelUnits] == crsUnits { + return nil + } + return errors.New("The unit system of the model and coordinate reference system are inconsistent") + } + return errors.New("Unable to check unit consistency, could not identify the coordinate reference system's units") +} + +func dataPairsfromTextBlock(sc *bufio.Scanner, nPairs int, colWidth int, valueWidth int) ([][2]float64, error) { + var stride int = valueWidth * 2 + pairs := [][2]float64{} +out: + for sc.Scan() { + line := sc.Text() + for s := 0; s < colWidth; { + if len(line) > s { + val1, err := strconv.ParseFloat(strings.TrimSpace(line[s:s+valueWidth]), 64) + if err != nil { + return pairs, err + } + val2, err := strconv.ParseFloat(strings.TrimSpace(line[s+valueWidth:s+stride]), 64) + if err != nil { + return pairs, err + } + pairs = append(pairs, [2]float64{val1, val2}) + if len(pairs) == nPairs { + break out + } + } else { + break + } + s += stride + } + } + return pairs, nil +} + +func getDataPairsfromTextBlock(nDataPairsLine string, sc *bufio.Scanner, colWidth int, valueWidth int) ([][2]float64, error) { + pairs := [][2]float64{} + for sc.Scan() { + line := sc.Text() + if strings.HasPrefix(line, nDataPairsLine) { + nPairs, err := strconv.Atoi(rightofEquals(line)) + if err != nil { + return pairs, err + } + pairs, err = dataPairsfromTextBlock(sc, nPairs, colWidth, valueWidth) + if err != nil { + return pairs, err + } + break + } + } + return pairs, nil +} + +// distance returns the distance along a straight line in euclidean space +func distance(p0, p1 [2]float64) float64 { + result := math.Sqrt(math.Pow((p1[0]-p0[0]), 2) + math.Pow((p1[1]-p0[1]), 2)) + return result +} + +// pointAtDistance returns a new point along a straight line in euclidean space +// at a specified distance +func pointAtDistance(p0, p1 [2]float64, delta float64) [2]float64 { + distanceRatio := delta / distance(p0, p1) + newX := (1-distanceRatio)*p0[0] + distanceRatio*p1[0] + newY := (1-distanceRatio)*p0[1] + distanceRatio*p1[1] + return [2]float64{newX, newY} +} + +// interpZ creates a new point a given distance along a line composed +// of many segments. +func interpXY(xyPairs [][2]float64, d float64) [2]float64 { + // newPoint is an x, y pair + var newPoint [2]float64 + lineSegments := len(xyPairs) - 1 + lineLength := 0.0 + +findLineSegment: + for i := 0; i < lineSegments; i++ { + p0, p1 := xyPairs[i], xyPairs[i+1] + lineLength += distance(p0, p1) + + switch { + case lineLength > d: + delta := distance(p0, p1) - (lineLength - d) + newPoint = pointAtDistance(p0, p1, delta) + break findLineSegment + + default: + continue + } + } + if d >= lineLength { + if d-lineLength <= 0.1 { + return xyPairs[len(xyPairs)-1] + } + fmt.Printf("The interpolated point has a station of %v while the xy line is %v long", d, lineLength) + } + return newPoint +} + +// attributeZ using station from cross-section line and gis coordinates +func attributeZ(xyPairs [][2]float64, mzPairs [][2]float64) []xyzPoint { + points := []xyzPoint{} + startingStation := mzPairs[0][0] + for _, mzPair := range mzPairs { + newPoint := interpXY(xyPairs, mzPair[0]-startingStation) + if newPoint[0] != 0 && newPoint[1] != 0 { + points = append(points, xyzPoint{newPoint[0], newPoint[1], mzPair[1]}) + } else { + fmt.Printf("Interpolated point has an xy value of (%v, %v). ", newPoint[0], newPoint[1]) + } + } + return points +} + +func getTransform(sourceCRS string, destinationCRS int) (gdal.CoordinateTransform, error) { + transform := gdal.CoordinateTransform{} + sourceSpRef := gdal.CreateSpatialReference(sourceCRS) + + destinationSpRef := gdal.CreateSpatialReference("") + if err := destinationSpRef.FromEPSG(destinationCRS); err != nil { + return transform, err + } + transform = gdal.CreateCoordinateTransform(sourceSpRef, destinationSpRef) + return transform, nil +} + +func flipXYLineString(xyLineString gdal.Geometry) gdal.Geometry { + yxLineString := gdal.Create(gdal.GT_LineString) + nPoints := xyLineString.PointCount() + for i := 0; i < nPoints; i++ { + x, y, _ := xyLineString.Point(i) + yxLineString.AddPoint2D(y, x) + } + xyLineString.Destroy() + return yxLineString +} + +func flipXYLineString25D(xyzLineString gdal.Geometry) gdal.Geometry { + yxzLineString := gdal.Create(gdal.GT_LineString25D) + nPoints := xyzLineString.PointCount() + for i := 0; i < nPoints; i++ { + x, y, z := xyzLineString.Point(i) + yxzLineString.AddPoint(y, x, z) + } + xyzLineString.Destroy() + return yxzLineString +} + +func flipXYLinearRing(xyLinearRing gdal.Geometry) gdal.Geometry { + yxLinearRing := gdal.Create(gdal.GT_LinearRing) + nPoints := xyLinearRing.PointCount() + for i := 0; i < nPoints; i++ { + x, y, _ := xyLinearRing.Point(i) + yxLinearRing.AddPoint2D(y, x) + } + xyLinearRing.Destroy() + return yxLinearRing +} + +func flipXYPoint(xyPoint gdal.Geometry) gdal.Geometry { + yxPoint := gdal.Create(gdal.GT_Point) + nPoints := xyPoint.PointCount() + for i := 0; i < nPoints; i++ { + x, y, _ := xyPoint.Point(i) + yxPoint.AddPoint2D(y, x) + } + xyPoint.Destroy() + return yxPoint +} + +func getRiverCenterline(sc *bufio.Scanner, transform gdal.CoordinateTransform) (VectorLayer, error) { + riverReach := strings.Split(rightofEquals(sc.Text()), ",") + layer := VectorLayer{FeatureName: fmt.Sprintf("%s, %s", strings.TrimSpace(riverReach[0]), strings.TrimSpace(riverReach[1]))} + + xyPairs, err := getDataPairsfromTextBlock("Reach XY=", sc, 64, 16) + if err != nil { + return layer, err + } + + xyLineString := gdal.Create(gdal.GT_LineString) + for _, pair := range xyPairs { + xyLineString.AddPoint2D(pair[0], pair[1]) + } + + xyLineString.Transform(transform) + // This is a temporary fix since the x and y values need to be flipped: + yxLineString := flipXYLineString(xyLineString) + + multiLineString := yxLineString.ForceToMultiLineString() + wkb, err := multiLineString.ToWKB() + if err != nil { + return layer, err + } + layer.Geometry = wkb + return layer, err +} + +func getXSBanks(sc *bufio.Scanner, transform gdal.CoordinateTransform, riverReachName string) (VectorLayer, []VectorLayer, error) { + bankLayers := []VectorLayer{} + + xsLayer, xyPairs, startingStation, err := getXS(sc, transform, riverReachName) + if err != nil { + return xsLayer, bankLayers, err + } + for sc.Scan() { + line := sc.Text() + if strings.HasPrefix(line, "Bank Sta=") { + bankLayers, err = getBanks(line, transform, xsLayer, xyPairs, startingStation) + if err != nil { + return xsLayer, bankLayers, err + } + break + } + } + return xsLayer, bankLayers, err +} + +func getXS(sc *bufio.Scanner, transform gdal.CoordinateTransform, riverReachName string) (VectorLayer, [][2]float64, float64, error) { + compData := strings.Split(rightofEquals(sc.Text()), ",") + layer := VectorLayer{FeatureName: strings.TrimSpace(compData[1]), Fields: map[string]interface{}{}} + layer.Fields["RiverReachName"] = riverReachName + + xyPairs, err := getDataPairsfromTextBlock("XS GIS Cut Line", sc, 64, 16) + if err != nil { + return layer, xyPairs, 0.0, err + } + + mzPairs, err := getDataPairsfromTextBlock("#Sta/Elev", sc, 80, 8) + if err != nil { + return layer, xyPairs, mzPairs[0][0], err + } + + xyzPoints := attributeZ(xyPairs, mzPairs) + + xyzLineString := gdal.Create(gdal.GT_LineString25D) + for _, point := range xyzPoints { + xyzLineString.AddPoint(point.x, point.y, point.z) + } + + xyzLineString.Transform(transform) + // This is a temporary fix since the x and y values need to be flipped + yxzLineString := flipXYLineString25D(xyzLineString) + + multiLineString := yxzLineString.ForceToMultiLineString() + wkb, err := multiLineString.ToWKB() + if err != nil { + return layer, xyPairs, mzPairs[0][0], err + } + layer.Geometry = wkb + return layer, xyPairs, mzPairs[0][0], err +} + +func getBanks(line string, transform gdal.CoordinateTransform, xsLayer VectorLayer, xyPairs [][2]float64, startingStation float64) ([]VectorLayer, error) { + layers := []VectorLayer{} + + bankStations := strings.Split(rightofEquals(line), ",") + for _, s := range bankStations { + layer := VectorLayer{FeatureName: strings.TrimSpace(s), Fields: map[string]interface{}{}} + layer.Fields["RiverReachName"] = xsLayer.Fields["RiverReachName"] + layer.Fields["xsName"] = xsLayer.FeatureName + bankStation, err := strconv.ParseFloat(s, 64) + if err != nil { + return layers, err + } + bankXY := interpXY(xyPairs, bankStation-startingStation) + 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 + yxPoint := flipXYPoint(xyPoint) + multiPoint := yxPoint.ForceToMultiPoint() + wkb, err := multiPoint.ToWKB() + if err != nil { + return layers, err + } + layer.Geometry = wkb + layers = append(layers, layer) + } + return layers, nil +} + +func getStorageArea(sc *bufio.Scanner, transform gdal.CoordinateTransform) (VectorLayer, error) { + layer := VectorLayer{FeatureName: strings.TrimSpace(strings.Split(rightofEquals(sc.Text()), ",")[0])} + + xyPairs, err := getDataPairsfromTextBlock("Storage Area Surface Line=", sc, 32, 16) + if err != nil { + return layer, err + } + + xyLinearRing := gdal.Create(gdal.GT_LinearRing) + for _, pair := range xyPairs { + xyLinearRing.AddPoint2D(pair[0], pair[1]) + } + + xyLinearRing.Transform(transform) + // This is a temporary fix since the x and y values need to be flipped: + yxLinearRing := flipXYLinearRing(xyLinearRing) + + yxPolygon := gdal.Create(gdal.GT_Polygon) + yxPolygon.AddGeometry(yxLinearRing) + yxMultiPolygon := yxPolygon.ForceToMultiPolygon() + wkb, err := yxMultiPolygon.ToWKB() + if err != nil { + return layer, err + } + layer.Geometry = wkb + return layer, err +} + +// GetGeospatialData ... +func GetGeospatialData(gd *GeoData, fs filestore.FileStore, geomFilePath string, sourceCRS string, destinationCRS int) error { + geomFileName := filepath.Base(geomFilePath) + f := Features{} + riverReachName := "" + + file, err := fs.GetObject(geomFilePath) + if err != nil { + return err + } + defer file.Close() + + sc := bufio.NewScanner(file) + + transform, err := getTransform(sourceCRS, destinationCRS) + if err != nil { + return err + } + + for sc.Scan() { + line := sc.Text() + switch { + case strings.HasPrefix(line, "River Reach="): + riverLayer, err := getRiverCenterline(sc, transform) + if err != nil { + return err + } + f.Rivers = append(f.Rivers, riverLayer) + riverReachName = riverLayer.FeatureName + + case strings.HasPrefix(line, "Storage Area="): + storageAreaLayer, err := getStorageArea(sc, transform) + if err != nil { + return err + } + f.StorageAreas = append(f.StorageAreas, storageAreaLayer) + + case strings.HasPrefix(line, "Type RM Length L Ch R = 1"): + xsLayer, bankLayers, err := getXSBanks(sc, transform, riverReachName) + if err != nil { + return err + } + f.XS = append(f.XS, xsLayer) + f.Banks = append(f.Banks, bankLayers...) + } + } + + gd.Features[geomFileName] = f + return nil +} diff --git a/tools/model.go b/tools/model.go index c897ebc..ac4d366 100644 --- a/tools/model.go +++ b/tools/model.go @@ -1,6 +1,7 @@ package tools import ( + "bufio" "errors" "fmt" "path/filepath" @@ -8,6 +9,7 @@ import ( "sync" "github.com/USACE/filestore" + "github.com/dewberry/gdal" ) type fileExtMatchers struct { @@ -21,6 +23,7 @@ type fileExtMatchers struct { SteadyRun *regexp.Regexp UnsteadyRun *regexp.Regexp AllFlowRun *regexp.Regexp + Projection *regexp.Regexp } var rasRE fileExtMatchers = fileExtMatchers{ // Maybe these ones are better? need a regex experts opinion @@ -34,13 +37,15 @@ var rasRE fileExtMatchers = fileExtMatchers{ // Maybe these ones are better? nee SteadyRun: regexp.MustCompile(".r[0-9][0-9]"), // `^\.r(0[1-9]|[1-9][0-9])$` UnsteadyRun: regexp.MustCompile(".x[0-9][0-9]"), // `^\.x(0[1-9]|[1-9][0-9])$` AllFlowRun: regexp.MustCompile(".[rx][0-9][0-9]"), // `^\.[rx](0[1-9]|[1-9][0-9])$` + Projection: regexp.MustCompile(".pr[oj]"), } // holder of multiple wait groups to help process files concurrency type rasWaitGroup struct { - Geom sync.WaitGroup - Plan sync.WaitGroup - Flow sync.WaitGroup + Geom sync.WaitGroup + Plan sync.WaitGroup + Flow sync.WaitGroup + Projection sync.WaitGroup } // Model is a general type should contain all necessary data for a model of any type. @@ -51,6 +56,7 @@ type Model struct { Files ModelFiles } +// ModelFiles ... type ModelFiles struct { InputFiles InputFiles OutputFiles OutputFiles @@ -100,6 +106,7 @@ type SupplementalFiles struct { ObservationalData interface{} // placeholder } +// RasModel ... type RasModel struct { FileStore filestore.FileStore ModelDirectory string @@ -110,10 +117,12 @@ type RasModel struct { FileList []string } +// IsAModel ... func (rm *RasModel) IsAModel() bool { return rm.isModel } +// IsGeospatial ... func (rm *RasModel) IsGeospatial() bool { if rm.Metadata.GeomFiles[0].FileExt != "" { return true @@ -121,14 +130,17 @@ func (rm *RasModel) IsGeospatial() bool { return false } +// ModelType ... func (rm *RasModel) ModelType() string { return rm.Type } +// ModelVersion ... func (rm *RasModel) ModelVersion() string { return rm.Version } +// Index ... func (rm *RasModel) Index() (Model, error) { if !rm.IsAModel() { return Model{}, errors.New("model is not valid") @@ -189,15 +201,31 @@ func (rm *RasModel) Index() (Model, error) { } // GeospatialData ... -func (rm *RasModel) GeospatialData() (GeoData, error) { - gd := GeoData{Features: make(map[string]Features), Georeference: DestinationCRS} +func (rm *RasModel) GeospatialData(destinationCRS int) (GeoData, error) { + gd := GeoData{} + + modelUnits := rm.Metadata.ProjFileContents.Units + + sourceCRS := rm.Metadata.Projection + + if sourceCRS == "" { + return gd, errors.New("Cannot extract geospatial data, no valid coordinate reference system") + } + + if err := checkUnitConsistency(modelUnits, sourceCRS); err != nil { + return gd, err + } + + gd.Features = make(map[string]Features) + gd.Georeference = destinationCRS for _, g := range rm.Metadata.GeomFiles { - if err := GetGeospatialData(&gd, rm.FileStore, g.Path); err != nil { + if err := GetGeospatialData(&gd, rm.FileStore, g.Path, sourceCRS, destinationCRS); err != nil { return gd, err } } return gd, nil + } func getModelFiles(rm *RasModel) error { @@ -209,14 +237,44 @@ func getModelFiles(rm *RasModel) error { } for _, file := range *files { - if filepath.Ext(file.Name) != ".prj" { - rm.FileList = append(rm.FileList, filepath.Join(file.Path, file.Name)) - } + rm.FileList = append(rm.FileList, filepath.Join(file.Path, file.Name)) } return nil } +// getProjection Reads a projection file. returns none to allow concurrency +func getProjection(rm *RasModel, fn string, wg *sync.WaitGroup, errChan chan error) { + + defer wg.Done() + + f, err := rm.FileStore.GetObject(fn) + if err != nil { + errChan <- err + return + } + defer f.Close() + + sc := bufio.NewScanner(f) + sc.Scan() + line := sc.Text() + + sourceSpRef := gdal.CreateSpatialReference(line) + if err := sourceSpRef.Validate(); err != nil { + rm.Metadata.Projection = "" + return + } + if rm.Metadata.Projection != "" { + errChan <- errors.New("Multiple projection files identified, cannot determine coordinate reference system") + return + } + + rm.Metadata.Projection = line + + return +} + +// NewRasModel ... func NewRasModel(key string, fs filestore.FileStore) (*RasModel, error) { rm := RasModel{ModelDirectory: filepath.Dir(key), FileStore: fs, Type: "RAS"} @@ -256,12 +314,19 @@ func NewRasModel(key string, fs filestore.FileStore) (*RasModel, error) { rasWG.Flow.Add(1) go getFlowData(&rm, fp, &rasWG.Flow, errChan) + case rasRE.Projection.MatchString(ext): + if filepath.Base(key) != filepath.Base(fp) { + rasWG.Projection.Add(1) + go getProjection(&rm, fp, &rasWG.Projection, errChan) + } + } } rasWG.Plan.Wait() rasWG.Geom.Wait() rasWG.Flow.Wait() + rasWG.Projection.Wait() if len(errChan) > 0 { fmt.Printf("Encountered %d errors\n", len(errChan)) @@ -269,13 +334,22 @@ func NewRasModel(key string, fs filestore.FileStore) (*RasModel, error) { } for _, p := range rm.Metadata.PlanFiles { - rm.Version += fmt.Sprintf("%s: %s, ", p.FileExt, p.ProgramVersion) + version := p.ProgramVersion + if version != "" { + rm.Version += fmt.Sprintf("%s: %s, ", p.FileExt, version) + } } for _, g := range rm.Metadata.GeomFiles { - rm.Version += fmt.Sprintf("%s: %s, ", g.FileExt, g.ProgramVersion) + version := g.ProgramVersion + if version != "" { + rm.Version += fmt.Sprintf("%s: %s, ", g.FileExt, version) + } } for _, f := range rm.Metadata.FlowFiles { - rm.Version += fmt.Sprintf("%s: %s, ", f.FileExt, f.ProgramVersion) + version := f.ProgramVersion + if version != "" { + rm.Version += fmt.Sprintf("%s: %s, ", f.FileExt, version) + } } if len(rm.Version) >= 2 { diff --git a/tools/project.go b/tools/project.go index b2567cc..c890d36 100644 --- a/tools/project.go +++ b/tools/project.go @@ -17,6 +17,7 @@ type ProjectMetadata struct { PlanFiles []PlanFileContents //`json:"Plan Data"` FlowFiles []FlowFileContents //`json:"Flow Data"` GeomFiles []GeomFileContents //`json:"Geometry Data"` + Projection string //`json:"Projection"` Notes string //`json:"Notes"` } @@ -52,7 +53,6 @@ func rmNewLineChar(s string) string { // verifyPrjPath identifies the .prj file within the passed model directory ... func verifyPrjPath(key string, rm *RasModel) error { - rasFiles := make([]string, 0) if filepath.Ext(key) != ".prj" { return fmt.Errorf("%s is not a .prj file", key) @@ -66,9 +66,7 @@ func verifyPrjPath(key string, rm *RasModel) error { return fmt.Errorf("%s is not a RAS Project file", key) } - rasFiles = append(rasFiles, key) rm.Metadata.ProjFilePath = key - rm.FileList = rasFiles return nil } diff --git a/tools/structures.go b/tools/structures.go new file mode 100644 index 0000000..12fe2d8 --- /dev/null +++ b/tools/structures.go @@ -0,0 +1,597 @@ +package tools + +import ( + "bufio" + "errors" + "math" + "strconv" + "strings" +) + +var conduitShapes map[int]string = map[int]string{ + 1: "Circular", + 2: "Box", + 3: "Pipe Arch", + 4: "Ellipse", + 5: "Arch", + 6: "Semi-Circle", + 7: "Low Arch", + 8: "High Arch", + 9: "Conspan Arch"} + +type hydraulicStructures struct { + River string `json:"River Name"` + Reach string `json:"Reach Name"` + NumXS int `json:"Num CrossSections"` + CulvertData culvertData `json:"Culvert Data"` + BridgeData bridgeData `json:"Bridge Data"` + WeirData weirData `json:"Inline Weir Data"` +} + +type culvertData struct { + NumCulverts int `json:"Num Culverts"` + Culverts []culverts `json:"Culverts"` +} + +type culverts struct { + Name string + Station float64 + Description string + DeckWidth float64 `json:"Deck Width"` + UpHighChord maxMinPairs `json:"Upstream High Chord"` + UpLowChord maxMinPairs `json:"Upstream Low Chord"` + DownHighChord maxMinPairs `json:"Downstream High Chord"` + DownLowChord maxMinPairs `json:"Downstream Low Chord"` + NumConduits int `json:"Num Culvert Conduits"` + Conduits []conduits `json:"Culvert Conduits"` +} + +type maxMinPairs struct { + Max float64 + Min float64 +} + +type conduits struct { + Name string + NumBarrels int `json:"Num Barrels"` + Shape string + Rise float64 + Span float64 + Length float64 + ManningsN float64 `json:"Mannings N"` +} + +type bridgeData struct { + NumBridges int `json:"Num Bridges"` + Bridges []bridges `json:"Bridges"` +} + +type bridges struct { + Name string + Station float64 + Description string + DeckWidth float64 `json:"Deck Width"` + UpHighChord maxMinPairs `json:"Upstream High Chord"` + UpLowChord maxMinPairs `json:"Upstream Low Chord"` + DownHighChord maxMinPairs `json:"Downstream High Chord"` + DownLowChord maxMinPairs `json:"Downstream Low Chord"` + NumPiers int `json:"Num Piers"` +} + +type weirData struct { + NumWeirs int `json:"Num Inline Weirs"` + Weirs []weirs `json:"Inline Weirs"` +} + +type weirs struct { + Name string + Station float64 + Description string + WeirWidth float64 `json:"Weir Width"` + WeirElev maxMinPairs `json:"Weir Elevations"` + NumGates int `json:"Num Gates"` + Gates []gates + NumConduits int `json:"Num Culvert Conduits"` + Conduits []conduits `json:"Culvert Conduits"` +} + +type gates struct { + Name string + Width float64 + Height float64 + NumOpenings int `json:"Num Openings"` +} + +func datafromTextBlock(hsSc *bufio.Scanner, i int, nLines int, nSkipLines int, colWidth int, valueWidth int, interval int) ([]float64, int, error) { + values := []float64{} + nSkipped := 0 + nProcessed := 0 + nvalues := 0 +out: + for hsSc.Scan() { + i++ + if nSkipped < nSkipLines { + nSkipped++ + continue + } + nProcessed++ + line := hsSc.Text() + for s := 0; s < colWidth; { + if len(line) > s { + sVal := strings.TrimSpace(line[s : s+valueWidth]) + if sVal != "" { + nvalues++ + if nvalues%interval == 0 { + val, err := strconv.ParseFloat(sVal, 64) + if err != nil { + return values, i, err + } + values = append(values, val) + } + } + s += valueWidth + } else { + if nLines == nProcessed { + break out + } + break + } + + } + if nLines == nProcessed { + break out + } + } + return values, i, nil +} + +func getMaxMinElev(hsSc *bufio.Scanner, i int, nLines int, nSkipLines int, colWidth int, valueWidth int, interval int) (maxMinPairs, int, error) { + pair := maxMinPairs{} + + elevations, i, err := datafromTextBlock(hsSc, i, nLines, nSkipLines, colWidth, valueWidth, interval) + + if err != nil { + return pair, i, err + } + + if len(elevations) == 0 { + return pair, i, nil + } + + maxElev, err := maxValue(elevations) + if err != nil { + return pair, i, err + } + + minElev, err := minValue(elevations) + if err != nil { + return pair, i, err + } + + pair = maxMinPairs{Max: maxElev, Min: minElev} + return pair, i, nil +} + +func numberofLines(nValues int, colWidth int, valueWidth int) int { + nLines := math.Ceil(float64(nValues) / (float64(colWidth) / float64(valueWidth))) + return int(nLines) +} + +func getHighLowChord(hsSc *bufio.Scanner, i int, nElevText string, colWidth int, valueWidth int) ([2]maxMinPairs, int, error) { + highLowPairs := [2]maxMinPairs{} + + nElev, err := strconv.Atoi(strings.TrimSpace(nElevText)) + if err != nil { + return highLowPairs, i, err + } + + nLines := numberofLines(nElev, 80, 8) + + highPair, i, err := getMaxMinElev(hsSc, i, nLines, nLines, 80, 8, 1) + if err != nil { + return highLowPairs, i, err + } + highLowPairs[0] = highPair + + lowPair, i, err := getMaxMinElev(hsSc, i, nLines, 0, 80, 8, 1) + if err != nil { + return highLowPairs, i, err + } + highLowPairs[1] = lowPair + + return highLowPairs, i, nil +} + +func stringtoFloat(s string) (float64, error) { + trimmed := strings.TrimSpace(s) + if trimmed != "" { + val, err := strconv.ParseFloat(trimmed, 64) + if err != nil { + return 0, err + } + return val, nil + } + return 0, nil +} + +func getConduits(line string, single bool) (conduits, error) { + lineData := strings.Split(rightofEquals(line), ",") + conduit := conduits{} + + if single { + conduit.NumBarrels = 1 + conduit.Name = strings.TrimSpace(lineData[13]) + + } else { + numbarrels, err := strconv.Atoi(strings.TrimSpace(lineData[11])) + if err != nil { + return conduit, err + } + conduit.NumBarrels = numbarrels + conduit.Name = strings.TrimSpace(lineData[12]) + } + + shapeID, err := strconv.Atoi(strings.TrimSpace(lineData[0])) + if err != nil { + return conduit, err + } + conduit.Shape = conduitShapes[shapeID] + + rise, err := stringtoFloat(lineData[1]) + if err != nil { + return conduit, err + } + conduit.Rise = rise + + span, err := stringtoFloat(lineData[2]) + if err != nil { + return conduit, err + } + conduit.Span = span + + length, err := stringtoFloat(lineData[3]) + if err != nil { + return conduit, err + } + conduit.Length = length + + mannings, err := stringtoFloat(lineData[4]) + if err != nil { + return conduit, err + } + conduit.ManningsN = mannings + + return conduit, nil +} + +func getCulvertData(hsSc *bufio.Scanner, i int, lineData []string) (culverts, int, error) { + culvert := culverts{} + + station, err := strconv.ParseFloat(strings.TrimSpace(lineData[1]), 64) + if err != nil { + return culvert, i, err + } + culvert.Station = station + + for hsSc.Scan() { + i++ + line := hsSc.Text() + switch { + case strings.HasPrefix(line, "BEGIN DESCRIPTION"): + var description string + description, i, err = getDescription(hsSc, i, "END DESCRIPTION:") + if err != nil { + return culvert, i, err + } + culvert.Description += description + + case strings.HasPrefix(line, "Node Name="): + culvert.Name = rightofEquals(line) + + case strings.HasPrefix(line, "Deck Dist"): + i++ + hsSc.Scan() + nextLineData := strings.Split(hsSc.Text(), ",") + deckWidth, err := strconv.ParseFloat(strings.TrimSpace(nextLineData[0]), 64) + if err != nil { + return culvert, i, err + } + culvert.DeckWidth = deckWidth + + var upHighLowPair [2]maxMinPairs + upHighLowPair, i, err = getHighLowChord(hsSc, i, nextLineData[4], 80, 8) + if err != nil { + return culvert, i, err + } + culvert.UpHighChord = upHighLowPair[0] + culvert.UpLowChord = upHighLowPair[1] + + var downHighLowPair [2]maxMinPairs + downHighLowPair, i, err = getHighLowChord(hsSc, i, nextLineData[5], 80, 8) + if err != nil { + return culvert, i, err + } + culvert.DownHighChord = downHighLowPair[0] + culvert.DownLowChord = downHighLowPair[1] + + case strings.HasPrefix(line, "Culvert="): + conduit, err := getConduits(line, true) + if err != nil { + return culvert, i, err + } + culvert.Conduits = append(culvert.Conduits, conduit) + culvert.NumConduits++ + + case strings.HasPrefix(line, "Multiple Barrel Culv="): + conduit, err := getConduits(line, false) + if err != nil { + return culvert, i, err + } + culvert.Conduits = append(culvert.Conduits, conduit) + culvert.NumConduits++ + + case strings.HasPrefix(line, "BC Design"): + return culvert, i, nil + + case strings.HasPrefix(line, "Type RM Length L Ch R ="): + return culvert, i, errors.New("Failed to terminate parsing of culvert at 'BC Design'") + + case strings.HasPrefix(line, "River Reach="): + return culvert, i, errors.New("Failed to terminate parsing of culvert at 'BC Design'") + } + } + return culvert, i, nil +} + +func getBridgeData(hsSc *bufio.Scanner, i int, lineData []string) (bridges, int, error) { + bridge := bridges{} + + station, err := strconv.ParseFloat(strings.TrimSpace(lineData[1]), 64) + if err != nil { + return bridge, i, err + } + bridge.Station = station + + for hsSc.Scan() { + i++ + line := hsSc.Text() + switch { + case strings.HasPrefix(line, "BEGIN DESCRIPTION"): + var description string + description, i, err = getDescription(hsSc, i, "END DESCRIPTION:") + if err != nil { + return bridge, i, err + } + bridge.Description += description + + case strings.HasPrefix(line, "Node Name="): + bridge.Name = rightofEquals(line) + + case strings.HasPrefix(line, "Deck Dist"): + i++ + hsSc.Scan() + nextLineData := strings.Split(hsSc.Text(), ",") + deckWidth, err := strconv.ParseFloat(strings.TrimSpace(nextLineData[0]), 64) + if err != nil { + return bridge, i, err + } + bridge.DeckWidth = deckWidth + + var upHighLowPair [2]maxMinPairs + upHighLowPair, i, err = getHighLowChord(hsSc, i, nextLineData[4], 80, 8) + if err != nil { + return bridge, i, err + } + bridge.UpHighChord = upHighLowPair[0] + bridge.UpLowChord = upHighLowPair[1] + + var downHighLowPair [2]maxMinPairs + downHighLowPair, i, err = getHighLowChord(hsSc, i, nextLineData[5], 80, 8) + if err != nil { + return bridge, i, err + } + bridge.DownHighChord = downHighLowPair[0] + bridge.DownLowChord = downHighLowPair[1] + + case strings.HasPrefix(line, "Pier Skew"): + bridge.NumPiers++ + + case strings.HasPrefix(line, "BC Design"): + return bridge, i, nil + + case strings.HasPrefix(line, "Type RM Length L Ch R ="): + return bridge, i, errors.New("Failed to terminate parsing of bridge at 'BC Design'") + + case strings.HasPrefix(line, "River Reach="): + return bridge, i, errors.New("Failed to terminate parsing of bridge at 'BC Design'") + } + } + return bridge, i, nil +} + +func getGates(nextLine string) (gates, error) { + gate := gates{} + + nextLineData := strings.Split(nextLine, ",") + + gate.Name = strings.TrimSpace(nextLineData[0]) + + width, err := stringtoFloat(nextLineData[1]) + if err != nil { + return gate, err + } + gate.Width = width + + height, err := stringtoFloat(nextLineData[2]) + if err != nil { + return gate, err + } + gate.Height = height + + numopenings, err := strconv.Atoi(strings.TrimSpace(nextLineData[13])) + if err != nil { + return gate, err + } + gate.NumOpenings = numopenings + + return gate, nil +} + +func getWeirData(rm *RasModel, fn string, i int) (weirs, error) { + weir := weirs{} + + f, err := rm.FileStore.GetObject(fn) + if err != nil { + return weir, err + } + defer f.Close() + + wSc := bufio.NewScanner(f) + + wi := 0 + for wSc.Scan() { + wi++ + if wi == i { + lineData := strings.Split(rightofEquals(wSc.Text()), ",") + station, err := strconv.ParseFloat(strings.TrimSpace(lineData[1]), 64) + if err != nil { + return weir, err + } + weir.Station = station + } else if wi > i { + line := wSc.Text() + switch { + case strings.HasPrefix(line, "BEGIN DESCRIPTION"): + description, _, err := getDescription(wSc, 0, "END DESCRIPTION:") + if err != nil { + return weir, err + } + weir.Description += description + + case strings.HasPrefix(line, "Node Name="): + weir.Name = rightofEquals(line) + + case strings.HasPrefix(line, "#Inline Weir SE="): + nElev, err := strconv.Atoi(strings.TrimSpace(rightofEquals(line))) + if err != nil { + return weir, err + } + nLines := numberofLines(nElev*2, 80, 8) + + elev, _, err := getMaxMinElev(wSc, 0, nLines, 0, 80, 8, 2) + if err != nil { + return weir, err + } + weir.WeirElev = elev + + case strings.HasPrefix(line, "IW Dist,WD"): + wSc.Scan() + nextLineData := strings.Split(wSc.Text(), ",") + weirWidth, err := strconv.ParseFloat(strings.TrimSpace(nextLineData[1]), 64) + if err != nil { + return weir, err + } + weir.WeirWidth = weirWidth + + case strings.HasPrefix(line, "IW Gate Name"): + wSc.Scan() + gate, err := getGates(wSc.Text()) + if err != nil { + return weir, err + } + weir.Gates = append(weir.Gates, gate) + weir.NumGates++ + + case strings.HasPrefix(line, "IW Culv="): + conduit, err := getConduits(line, false) + if err != nil { + return weir, err + } + weir.Conduits = append(weir.Conduits, conduit) + weir.NumConduits++ + + case strings.HasPrefix(line, "Type RM Length L Ch R ="): + return weir, nil + + case strings.HasPrefix(line, "River Reach="): + return weir, nil + } + } + } + return weir, nil +} + +func getHydraulicStructureData(rm *RasModel, fn string, idx int) (hydraulicStructures, error) { + structures := hydraulicStructures{} + bData := bridgeData{} + cData := culvertData{} + wData := weirData{} + + f, err := rm.FileStore.GetObject(fn) + if err != nil { + return structures, err + } + defer f.Close() + + hsSc := bufio.NewScanner(f) + + i := 0 + for hsSc.Scan() { + i++ + if i == idx { + riverReach := strings.Split(rightofEquals(hsSc.Text()), ",") + structures.River = strings.TrimSpace(riverReach[0]) + structures.Reach = strings.TrimSpace(riverReach[1]) + } else if i > idx { + line := hsSc.Text() + if strings.HasPrefix(line, "Type RM Length L Ch R =") { + data := strings.Split(rightofEquals(line), ",") + structureType, err := strconv.Atoi(strings.TrimSpace(data[0])) + if err != nil { + return structures, err + } + switch structureType { + case 1: + structures.NumXS++ + + case 2: + var culvert culverts + culvert, i, err = getCulvertData(hsSc, i, data) + if err != nil { + return structures, err + } + cData.Culverts = append(cData.Culverts, culvert) + cData.NumCulverts++ + + case 3: + var bridge bridges + bridge, i, err = getBridgeData(hsSc, i, data) + if err != nil { + return structures, err + } + bData.Bridges = append(bData.Bridges, bridge) + bData.NumBridges++ + + case 5: + weir, err := getWeirData(rm, fn, i) + if err != nil { + return structures, err + } + wData.Weirs = append(wData.Weirs, weir) + wData.NumWeirs++ + } + } + if strings.HasPrefix(line, "River Reach=") { + structures.CulvertData = cData + structures.BridgeData = bData + structures.WeirData = wData + return structures, nil + } + } + } + structures.CulvertData = cData + structures.BridgeData = bData + structures.WeirData = wData + + return structures, nil +} diff --git a/tools/utils.go b/tools/utils.go new file mode 100644 index 0000000..f404458 --- /dev/null +++ b/tools/utils.go @@ -0,0 +1,61 @@ +package tools + +import ( + "bufio" + "errors" + "strings" +) + +func maxValue(values []float64) (float64, error) { + if len(values) == 0 { + return 0.0, errors.New("Cannot detect a maximum value in an empty slice") + } + + max := values[0] + for _, v := range values { + if v > max { + max = v + } + } + + return max, nil +} + +func minValue(values []float64) (float64, error) { + if len(values) == 0 { + return 0.0, errors.New("Cannot detect a minimum value in an empty slice") + } + + min := values[0] + for _, v := range values { + if v < min { + min = v + } + } + + return min, nil +} + +func rightofEquals(line string) string { + return strings.TrimSpace(strings.Split(line, "=")[1]) +} + +func getDescription(sc *bufio.Scanner, idx int, endLine string) (string, int, error) { + description := "" + nLines := 0 + for sc.Scan() { + idx++ + line := sc.Text() + if strings.HasPrefix(line, endLine) { + return description, idx, nil + } + if line != "" { + if nLines > 0 { + description += "\n" + } + description += line + nLines++ + } + } + return description, idx, nil +}