Geospatial

Consequences.sln

@@ -14,7 +14,9 @@ Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "Solution Items", "Solution
 EndProject
 Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "ConsequencesTest", "ConsequencesTest\ConsequencesTest.csproj", "{95BCC49B-7780-41E9-8365-C51B5E1B3D5E}"
 EndProject
-Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "ScratchPaper", "ScratchPaper\ScratchPaper.csproj", "{31071E1B-DA08-40CF-8F16-95982E85B45D}"
+Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "ScratchPaper", "ScratchPaper\ScratchPaper.csproj", "{31071E1B-DA08-40CF-8F16-95982E85B45D}"
+EndProject
+Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "Geospatial", "Geospatial\Geospatial.csproj", "{CB7204A4-C01A-4C1C-9BEB-C82F95EEA214}"
 EndProject
 Global
 	GlobalSection(SolutionConfigurationPlatforms) = preSolution
@@ -34,6 +36,10 @@ Global
 		{31071E1B-DA08-40CF-8F16-95982E85B45D}.Debug|Any CPU.Build.0 = Debug|Any CPU
 		{31071E1B-DA08-40CF-8F16-95982E85B45D}.Release|Any CPU.ActiveCfg = Release|Any CPU
 		{31071E1B-DA08-40CF-8F16-95982E85B45D}.Release|Any CPU.Build.0 = Release|Any CPU
+		{CB7204A4-C01A-4C1C-9BEB-C82F95EEA214}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
+		{CB7204A4-C01A-4C1C-9BEB-C82F95EEA214}.Debug|Any CPU.Build.0 = Debug|Any CPU
+		{CB7204A4-C01A-4C1C-9BEB-C82F95EEA214}.Release|Any CPU.ActiveCfg = Release|Any CPU
+		{CB7204A4-C01A-4C1C-9BEB-C82F95EEA214}.Release|Any CPU.Build.0 = Release|Any CPU
 	EndGlobalSection
 	GlobalSection(SolutionProperties) = preSolution
 		HideSolutionNode = FALSE

Consequences/Consequences/IFileStreamingProcessor.cs

@@ -0,0 +1,7 @@
+using USACE.HEC.Geography;
+
+namespace USACE.HEC.Consequences;
+public interface IFileStreamingProcessor
+{
+  public void Process(string filePath, Action<IConsequencesReceptor> consequenceReceptorProcess);
+}

Consequences/Consequences/Structure.cs

@@ -103,4 +103,28 @@ public Result Compute(IHazard hazard)
 
     return new Result([.. resultItems]);
   }
+
+  public Result ToResult()
+  {
+    List<ResultItem> resultItems = [];
+    resultItems.Add(new ResultItem { ResultName = "fd_id", Result = Name });
+    resultItems.Add(new ResultItem { ResultName = "occtype", Result = Occtype });
+    resultItems.Add(new ResultItem { ResultName = "st_damcat", Result = DamCat });
+    resultItems.Add(new ResultItem { ResultName = "bldgtype", Result = ConstructionType });
+    resultItems.Add(new ResultItem { ResultName = "found_type", Result = FoundationType });
+    resultItems.Add(new ResultItem { ResultName = "cbfips", Result = CBFips });
+    resultItems.Add(new ResultItem { ResultName = "pop2amu65", Result = Popu65night });
+    resultItems.Add(new ResultItem { ResultName = "pop2amo65", Result = Popo65night });
+    resultItems.Add(new ResultItem { ResultName = "pop2pmu65", Result = Popu65day });
+    resultItems.Add(new ResultItem { ResultName = "pop2pmo65", Result = Popo65day });
+    resultItems.Add(new ResultItem { ResultName = "num_story", Result = NumStories });
+    resultItems.Add(new ResultItem { ResultName = "found_ht", Result = FoundHt });
+    resultItems.Add(new ResultItem { ResultName = "val_struct", Result = StructVal });
+    resultItems.Add(new ResultItem { ResultName = "val_cont", Result = ContVal });
+    resultItems.Add(new ResultItem { ResultName = "firmzone", Result = FirmZone });
+    resultItems.Add(new ResultItem { ResultName = "x", Result = X });
+    resultItems.Add(new ResultItem { ResultName = "y", Result = Y });
+    resultItems.Add(new ResultItem { ResultName = "ground_elv", Result = GroundElevation });
+    return new Result([.. resultItems]);
+  }
 }

Geospatial/Geospatial.csproj

@@ -0,0 +1,32 @@
+<Project Sdk="Microsoft.NET.Sdk">
+
+  <PropertyGroup>
+    <TargetFramework>net8.0</TargetFramework>
+    <ImplicitUsings>enable</ImplicitUsings>
+    <Nullable>enable</Nullable>
+    <OutputType>Library</OutputType>
+  </PropertyGroup>
+
+  <ItemGroup>
+	  <PackageReference Include="Geospatial.GDALAssist" Version="0.1.0-Beta" />
+  </ItemGroup>
+
+  <!--<ItemGroup>
+    <PackageReference Include="Geospatial.GDALAssist" Version="1.0.845-beta" />
+  </ItemGroup>-->
+
+
+
+  <ItemGroup>
+    <ProjectReference Include="..\Consequences\Consequences.csproj" />
+  </ItemGroup>
+
+
+
+  <!--<ItemGroup>
+    <Reference Include="Geospatial.GDALAssist">
+      <HintPath>..\..\..\Users\HEC\Downloads\GDAL\Geospatial.GDALAssist.dll</HintPath>
+    </Reference>
+  </ItemGroup>-->
+
+</Project>

Geospatial/SpatialStructureProcessor.cs

@@ -0,0 +1,45 @@
+using System.Reflection;
+using System.Text.Json.Serialization;
+using OSGeo.OGR;
+using USACE.HEC.Consequences;
+
+namespace Geospatial;
+
+// process an OGR driver into a stream of structures and apply a process to the structure
+public class SpatialStructureProcessor : IFileStreamingProcessor
+{
+  public void Process(string filePath, Action<IConsequencesReceptor> consequenceReceptorProcess)
+  {
+    GDALAssist.GDALSetup.InitializeMultiplatform();
+
+    using DataSource ds = Ogr.Open(filePath, 0) ?? throw new Exception("Failed to create datasource.");
+    Layer layer = ds.GetLayerByIndex(0) ?? throw new Exception("Failed to create layer.");
+
+    Feature structure;
+    while ((structure = layer.GetNextFeature()) != null)
+    {
+      PropertyInfo[] properties = typeof(Structure).GetProperties();
+      Structure s = new();
+
+      foreach (PropertyInfo property in properties)
+      {
+        JsonPropertyNameAttribute? jsonPropertyAttribute = property.GetCustomAttribute<JsonPropertyNameAttribute>();
+
+        // check if the property has a JsonPropertyName (jsonPropertyAttribute not null) 
+        // we know that all strucuture properties have one but the compiler does not?
+        string jsonPropertyName = jsonPropertyAttribute != null ? jsonPropertyAttribute.Name : property.Name;
+
+        if (property.PropertyType == typeof(int))
+          property.SetValue(s, structure.GetFieldAsInteger(jsonPropertyName));
+        else if (property.PropertyType == typeof(double)) 
+          property.SetValue(s, structure.GetFieldAsDouble(jsonPropertyName));
+        else if (property.PropertyType == typeof(float))
+          property.SetValue(s, (float)structure.GetFieldAsDouble(jsonPropertyName));
+        else // field is a string
+          property.SetValue(s, structure.GetFieldAsString(jsonPropertyName));
+      }
+
+      consequenceReceptorProcess(s);
+    }
+  }
+}

Geospatial/SpatialWriter.cs

@@ -0,0 +1,112 @@
+using OSGeo.OGR;
+using OSGeo.OSR;
+using USACE.HEC.Results;
+
+namespace Geospatial;
+
+
+public class SpatialWriter : IResultsWriter
+{
+  private string _outputPath;
+  private Layer? _layer;
+  private DataSource? _dataSource;
+  private SpatialReference? _srs;
+  private SpatialReference? _dst;
+  private bool _headersWritten;
+  public delegate void FieldTypeDelegate(Feature layer, string fieldName, object value);
+  private FieldTypeDelegate? _fieldTypeDelegate;
+
+  public SpatialWriter(string outputPath, string driverName, int projection, FieldTypeDelegate fieldTypeDelegate)
+  {
+    GDALAssist.GDALSetup.InitializeMultiplatform();
+    _outputPath = outputPath;
+    _dataSource = Ogr.GetDriverByName(driverName).CreateDataSource(_outputPath, null);
+    if (_dataSource == null) throw new Exception("Failed to create data source.");
+    _srs = new SpatialReference("");
+    _srs.SetWellKnownGeogCS("WGS84");
+    if (_srs == null) throw new Exception("Failed to create SpatialReference.");
+    _dst = new SpatialReference("");
+    _dst.ImportFromEPSG(projection);
+    if (_dst == null) throw new Exception("Failed to create SpatialReference.");
+
+    _layer = _dataSource.CreateLayer("layer_name", _dst, wkbGeometryType.wkbPoint, null);
+    if (_layer == null) throw new Exception("Failed to create layer.");
+
+    _headersWritten = false;
+    _fieldTypeDelegate = fieldTypeDelegate;
+  }
+
+  public void Write(Result res)
+  {
+    if (_layer == null || _fieldTypeDelegate == null)
+    {
+      return;
+    }
+    if (!_headersWritten)
+    {
+      InitializeFields(_layer, res);
+      _headersWritten= true;
+    }
+
+    using Feature feature = new Feature(_layer.GetLayerDefn());
+    using Geometry geometry = new Geometry(wkbGeometryType.wkbPoint);
+    double x = (double)res.Fetch("x").Result;
+    double y = (double)res.Fetch("y").Result;
+    geometry.AddPoint(y, x, 0); 
+    // transform the coordinates according to the specified projection
+    CoordinateTransformation ct = new CoordinateTransformation(_srs, _dst);
+    geometry.Transform(ct);
+    feature.SetGeometry(geometry);
+
+    for (int i = 0; i < _layer.GetLayerDefn().GetFieldCount(); i++)
+    {
+      string fieldName = _layer.GetLayerDefn().GetFieldDefn(i).GetName();
+      object val = res.Fetch(fieldName).Result;
+      // assign value to field according to the result's field type mappings defined in _fieldTypeDelegate
+      _fieldTypeDelegate(feature, fieldName, val);
+    }
+    _layer.CreateFeature(feature);
+  }
+
+
+  // create fields of the appropriate types
+  private static void InitializeFields(Layer layer, Result res)
+  {
+    foreach (ResultItem item in res.ResultItems)
+    {
+      switch (item.Result)
+      {
+        case int _:
+          layer.CreateField(new FieldDefn(item.ResultName, FieldType.OFTInteger), 1);
+          break;
+        case long _:
+          layer.CreateField(new FieldDefn(item.ResultName, FieldType.OFTInteger64), 1);
+          break;
+        case double _ or float _:
+          layer.CreateField(new FieldDefn(item.ResultName, FieldType.OFTReal), 1);
+          break;
+        case string _:
+          layer.CreateField(new FieldDefn(item.ResultName, FieldType.OFTString), 1);
+          break;
+        case DateOnly _:
+          layer.CreateField(new FieldDefn(item.ResultName, FieldType.OFTDate), 1);
+          break;
+        case TimeOnly _:
+          layer.CreateField(new FieldDefn(item.ResultName, FieldType.OFTTime), 1);
+          break;
+        case DateTime _:
+          layer.CreateField(new FieldDefn(item.ResultName, FieldType.OFTDateTime), 1);
+          break;
+        default:
+          // for case of a null string
+          layer.CreateField(new FieldDefn(item.ResultName, FieldType.OFTString), 1);
+          break;
+      }
+    }
+  }
+
+  public void Dispose()
+  {
+
+  }
+}

Geospatial/Utilities.cs

@@ -0,0 +1,23 @@
+using OSGeo.OGR;
+using OSGeo.OSR;
+
+namespace Geospatial;
+public class Utilities
+{
+  public static void StructureFieldTypes(Feature feature, string fieldName, object value)
+  {
+    switch (fieldName)
+    {
+      case "fd_id" or "num_story" or "pop2pmo65" or "pop2pmu65" or "pop2amo65" or "pop2amu65":
+        feature.SetField(fieldName, (int)value);
+        break;
+      case "x" or "y" or "ground_elv" or "val_struct" or "val_cont" or "found_ht":
+        feature.SetField(fieldName, (double)value);
+        break;
+      case "st_damcat" or "cbfips" or "occtype" or "found_type" or "firmzone" or "bldgtype":
+        feature.SetField(fieldName, (string)value);
+        break;
+    }
+  }
+}
+

ScratchPaper/Program.cs

@@ -1,33 +1,51 @@
 using USACE.HEC.Consequences;
 using USACE.HEC.Geography;
+using Geospatial;
+using USACE.HEC.Results;
+using USACE.HEC.Hazards;
+using OSGeo.OGR;
 
 internal class Program
 {
-  private static async Task Main(string[] args)
+  private async static Task Main(string[] args)
   {
-    // city block in Sunset District, SF
-    Location upperLeft1 = new Location { X = -122.475275, Y = 37.752394 };
-    Location lowerRight1 = new Location { X = -122.473523, Y = 37.750642 };
+    // city blocks in Sunset District, SF
+    Location upperLeft1 = new Location { X = -122.48, Y = 37.76 };
+    Location lowerRight1 = new Location { X = -122.47, Y = 37.75 };
     BoundingBox boundingBox1 = new BoundingBox(upperLeft1, lowerRight1);
 
 
-    Location upperLeft2 = new Location { X = -122.5, Y = 37.8 };
-    Location lowerRight2 = new Location { X = -122, Y = 37.3 };
+    Location upperLeft2 = new Location { X = -121.74, Y = 38.58 };
+    Location lowerRight2 = new Location { X = -121.70, Y = 38.54 };
     BoundingBox boundingBox2 = new BoundingBox(upperLeft2, lowerRight2);
 
-    IBBoxStreamingProcessor sp = new NSIStreamingProcessor();
+    NSIStreamingProcessor sp = new NSIStreamingProcessor();
+    string filePath = @"C:\repos\consequences\ScratchPaper\generated";
 
+    //using SpatialWriter c = new SpatialWriter(filePath, "ESRI Shapefile", 3310, Utilities.StructureFieldTypes);
     int count = 0;
-    var watch = System.Diagnostics.Stopwatch.StartNew();
+
     await sp.Process(boundingBox2, (IConsequencesReceptor s) => {
       //Console.WriteLine(((Structure)s).Name);
+      Result res = ((Structure)s).ToResult();
+      //c.Write(res);
       count++;
     });
-    watch.Stop();
-    var elapsedMs = watch.ElapsedMilliseconds;
-
     Console.WriteLine(count);
-    Console.WriteLine("Time elapsed: " + elapsedMs.ToString() + " milliseconds");
-    Console.Read(); 
+    Console.Read();
+
+  }
+
+  public static void Read()
+  {
+    SpatialStructureProcessor reader = new SpatialStructureProcessor();
+    string path = @"C:\Data\Muncie_WS6_Solution_PART2\Muncie_WS6_Part1_Solution_PART2\Muncie_WS6_Part1_Solution\Structure Inventories\Existing_BaseSI\BaseMuncieStructsFinal.shp";
+    int count = 0;
+    reader.Process(path, (IConsequencesReceptor s) => {
+      Console.WriteLine($"Structure {count}:");
+      Console.WriteLine($"  fd_id: {((Structure)s).Name}");
+      Console.WriteLine($"  cbfips: {((Structure)s).CBFips}");
+      count++;
+    });
   }
 }

ScratchPaper/ScratchPaper.csproj

@@ -9,6 +9,7 @@
 
   <ItemGroup>
     <ProjectReference Include="..\Consequences\Consequences.csproj" />
+    <ProjectReference Include="..\Geospatial\Geospatial.csproj" />
   </ItemGroup>
 
 </Project>