From 8dc2a071095aaf62df2e1473ba0361f69fb02459 Mon Sep 17 00:00:00 2001 From: Shane Date: Fri, 12 Feb 2021 12:10:33 -0500 Subject: [PATCH] Update isgeospatial and geospatialdata for RAS version less than 4.0 --- tools/geospatial.go | 31 ++++++++++++++++++++------- tools/model.go | 51 ++++++++++++++++++++++++++++++--------------- 2 files changed, 58 insertions(+), 24 deletions(-) diff --git a/tools/geospatial.go b/tools/geospatial.go index 7bdbf9c..7302f20 100644 --- a/tools/geospatial.go +++ b/tools/geospatial.go @@ -4,6 +4,7 @@ import ( "bufio" "errors" "fmt" + "log" "math" "path/filepath" "regexp" @@ -271,6 +272,7 @@ func getXSBanks(sc *bufio.Scanner, transform gdal.CoordinateTransform, riverReac if err != nil { return xsLayer, bankLayers, err } + log.Println("Extracted cross-section") if xsLayer.Fields["CutLineProfileMatch"].(bool) { for sc.Scan() { line := sc.Text() @@ -283,6 +285,7 @@ func getXSBanks(sc *bufio.Scanner, transform gdal.CoordinateTransform, riverReac } } } + log.Println("Extracted banks") return xsLayer, bankLayers, err } @@ -304,6 +307,12 @@ func getXS(sc *bufio.Scanner, transform gdal.CoordinateTransform, riverReachName if err != nil { return layer, xyPairs, 0.0, err } + + if len(xyPairs) < 2 { + err = errors.New("the cross-section cutline could not be extracted, check that the geometry file contains cutlines") + return layer, xyPairs, 0.0, err + } + xyzLineString := gdal.Create(gdal.GT_LineString25D) for _, pair := range xyPairs { xyzLineString.AddPoint(pair[0], pair[1], 0.0) @@ -314,15 +323,17 @@ func getXS(sc *bufio.Scanner, transform gdal.CoordinateTransform, riverReachName if err != nil { return layer, xyPairs, mzPairs[0][0], err } - lenProfile := mzPairs[len(mzPairs)-1][0] - mzPairs[0][0] - if math.Abs(lenProfile-lenCutLine) <= 0.1 { - xyzPoints := attributeZ(xyPairs, mzPairs) - xyzLineString = gdal.Create(gdal.GT_LineString25D) - for _, point := range xyzPoints { - xyzLineString.AddPoint(point.x, point.y, point.z) + if len(mzPairs) >= 2 { + lenProfile := mzPairs[len(mzPairs)-1][0] - mzPairs[0][0] + if math.Abs(lenProfile-lenCutLine) <= 0.1 { + xyzPoints := attributeZ(xyPairs, mzPairs) + xyzLineString = gdal.Create(gdal.GT_LineString25D) + for _, point := range xyzPoints { + xyzLineString.AddPoint(point.x, point.y, point.z) + } + layer.Fields["CutLineProfileMatch"] = true } - layer.Fields["CutLineProfileMatch"] = true } xyzLineString.Transform(transform) @@ -400,6 +411,7 @@ func GetGeospatialData(gd *GeoData, fs filestore.FileStore, geomFilePath string, geomFileName := filepath.Base(geomFilePath) f := Features{} riverReachName := "" + log.Println("Extracting geospatial data from:", geomFilePath) file, err := fs.GetObject(geomFilePath) if err != nil { @@ -416,6 +428,7 @@ func GetGeospatialData(gd *GeoData, fs filestore.FileStore, geomFilePath string, for sc.Scan() { line := sc.Text() + switch { case strings.HasPrefix(line, "River Reach="): riverLayer, err := getRiverCenterline(sc, transform) @@ -424,6 +437,7 @@ func GetGeospatialData(gd *GeoData, fs filestore.FileStore, geomFilePath string, } f.Rivers = append(f.Rivers, riverLayer) riverReachName = riverLayer.FeatureName + log.Println("Extracted river centerline") case strings.HasPrefix(line, "Storage Area="): storageAreaLayer, err := getStorageArea(sc, transform) @@ -431,6 +445,7 @@ func GetGeospatialData(gd *GeoData, fs filestore.FileStore, geomFilePath string, return err } f.StorageAreas = append(f.StorageAreas, storageAreaLayer) + log.Println("Extracted storage area") case strings.HasPrefix(line, "Type RM Length L Ch R = 1"): xsLayer, bankLayers, err := getXSBanks(sc, transform, riverReachName) @@ -439,6 +454,8 @@ func GetGeospatialData(gd *GeoData, fs filestore.FileStore, geomFilePath string, } f.XS = append(f.XS, xsLayer) f.Banks = append(f.Banks, bankLayers...) + log.Println("Extracted banks and cross-sections") + } } diff --git a/tools/model.go b/tools/model.go index c53e506..ac28137 100644 --- a/tools/model.go +++ b/tools/model.go @@ -4,8 +4,11 @@ import ( "bufio" "errors" "fmt" + "log" "path/filepath" "regexp" + "strconv" + "strings" "sync" "github.com/USACE/filestore" @@ -196,9 +199,24 @@ func (rm *RasModel) Index() Model { // IsGeospatial ... func (rm *RasModel) IsGeospatial() bool { if rm.Metadata.Projection == "" { - fmt.Println("no valid coordinate reference system") + log.Println("no valid coordinate reference system") return false } + modelVersions := strings.Split(rm.Version, ",") + for _, version := range modelVersions { + if strings.Contains(version, ".g") { + geomVersion := strings.TrimSpace(strings.Split(version, ":")[1]) + v, err := strconv.ParseFloat(geomVersion, 64) + if err != nil { + log.Println("could not convert the geometry version to a float") + return false + } + if v < 4 { + log.Printf("geometry file version: %f is not geospatial", v) + return false + } + } + } return true } @@ -206,28 +224,27 @@ func (rm *RasModel) IsGeospatial() bool { // GeospatialData ... func (rm *RasModel) GeospatialData(destinationCRS int) (GeoData, error) { gd := GeoData{} + if rm.IsGeospatial() { + modelUnits := rm.Metadata.ProjFileContents.Units - modelUnits := rm.Metadata.ProjFileContents.Units - - sourceCRS := rm.Metadata.Projection - - if sourceCRS == "" { - return gd, errors.New("Cannot extract geospatial data, no valid coordinate reference system") - } + sourceCRS := rm.Metadata.Projection - if err := checkUnitConsistency(modelUnits, sourceCRS); err != nil { - return gd, err - } + if err := checkUnitConsistency(modelUnits, sourceCRS); err != nil { + return gd, err + } - gd.Features = make(map[string]Features) - gd.Georeference = destinationCRS + gd.Features = make(map[string]Features) + gd.Georeference = destinationCRS - for _, g := range rm.Metadata.GeomFiles { - if err := GetGeospatialData(&gd, rm.FileStore, g.Path, sourceCRS, destinationCRS); err != nil { - return gd, err + for _, g := range rm.Metadata.GeomFiles { + if err := GetGeospatialData(&gd, rm.FileStore, g.Path, sourceCRS, destinationCRS); err != nil { + return gd, err + } } + return gd, nil } - return gd, nil + err := errors.New("the model is not geospatial") + return gd, err }