Geospatial (#30)

* test opening a gpkg file

* base functionality for SpatialWriter

* initial OGR structure reader implementation

* renamed to SpatialStructureProcessor

* messing with different gdal build

* improvements to SpatialWriter

* added structure writing capabilities

* refactored SpatialProcessor to be more concise

* added constructor to SpatialProcessor. created InitializeGDAL function

* made SpatialProcessor a generic

* integrated reflection and generics into ToResult

---------

Co-authored-by: Brennan Beam <[email protected]>

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/IStreamingProcessor.cs

@@ -0,0 +1,7 @@
+using USACE.HEC.Geography;
+
+namespace USACE.HEC.Consequences;
+public interface IStreamingProcessor
+{
+  public void Process<T>(Action<IConsequencesReceptor> consequenceReceptorProcess) where T : IConsequencesReceptor, new();
+}

Consequences/Results/Utilities.cs

@@ -0,0 +1,24 @@
+using System.Reflection;
+using System.Text.Json.Serialization;
+using USACE.HEC.Consequences;
+
+namespace USACE.HEC.Results;
+public class Utilities
+{
+  public static Result ConsequenceReceptorToResult<T>(IConsequencesReceptor cr) where T: IConsequencesReceptor
+  { 
+    List<ResultItem> resultItems = [];
+
+    PropertyInfo[] properties = typeof(T).GetProperties();
+    foreach (PropertyInfo property in properties)
+    {
+      ResultItem item; 
+      JsonPropertyNameAttribute jsonPropertyAttribute = property.GetCustomAttribute<JsonPropertyNameAttribute>();
+      item.ResultName = jsonPropertyAttribute != null ? jsonPropertyAttribute.Name : property.Name;
+      item.Result = property.GetValue(cr);
+      resultItems.Add(item);
+    }
+
+    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/SpatialProcessor.cs

@@ -0,0 +1,50 @@
+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 IConsequenceReceptors and apply a consequenceReceptorProcess to each
+public class SpatialProcessor : IStreamingProcessor
+{
+  private DataSource _dataSource;
+  private Layer _layer;
+  public SpatialProcessor(string filePath)
+  {
+    _dataSource = Ogr.Open(filePath, 0) ?? throw new Exception("Failed to create datasource.");
+    _layer = _dataSource.GetLayerByIndex(0) ?? throw new Exception("Failed to create layer.");
+  }
+  public void Process<T>(Action<IConsequencesReceptor> consequenceReceptorProcess) where T : IConsequencesReceptor, new()
+  {
+    Feature feature;
+    while ((feature = _layer.GetNextFeature()) != null)
+    {
+      PropertyInfo[] properties = typeof(T).GetProperties();
+      // T MUST be an IConsequencesReceptor with a parameterless constructor, eg. a structure
+      T consequenceReceptor = new();
+
+      foreach (PropertyInfo property in properties)
+      {
+        // IConsequenceReceptors' properties must have the JsonPropertyName tag
+        // JsonPropertyNames must match the corresponding field names in the driver for a given property
+        JsonPropertyNameAttribute? jsonPropertyAttribute = property.GetCustomAttribute<JsonPropertyNameAttribute>();
+        string fieldName = jsonPropertyAttribute != null ? jsonPropertyAttribute.Name : property.Name;
+
+        // read values from the driver into their corresponding properties in the IConsequencesReceptor
+        if (property.PropertyType == typeof(int))
+          property.SetValue(consequenceReceptor, feature.GetFieldAsInteger(fieldName));
+        else if (property.PropertyType == typeof(long))
+          property.SetValue(consequenceReceptor, feature.GetFieldAsInteger64(fieldName));
+        else if (property.PropertyType == typeof(double)) 
+          property.SetValue(consequenceReceptor, feature.GetFieldAsDouble(fieldName));
+        else if (property.PropertyType == typeof(float))
+          property.SetValue(consequenceReceptor, (float)feature.GetFieldAsDouble(fieldName));
+        else if (property.PropertyType == typeof(string))
+          property.SetValue(consequenceReceptor, feature.GetFieldAsString(fieldName));
+      }
+
+      consequenceReceptorProcess(consequenceReceptor);
+    }
+  }
+}

Geospatial/SpatialWriter.cs

@@ -0,0 +1,107 @@
+using OSGeo.OGR;
+using OSGeo.OSR;
+using USACE.HEC.Results;
+
+namespace Geospatial;
+
+
+public class SpatialWriter : IResultsWriter
+{
+  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)
+  {
+    _dataSource = Ogr.GetDriverByName(driverName).CreateDataSource(outputPath, null) ?? throw new Exception("Failed to create data source.");
+    _srs = new SpatialReference("");
+    _srs.ImportFromEPSG(4326); // 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) ?? 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,28 @@
+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;
+    }
+  }
+
+  public static void InitializeGDAL()
+  {
+    GDALAssist.GDALSetup.InitializeMultiplatform();
+  }
+}
+

ScratchPaper/Program.cs

@@ -1,33 +1,54 @@
 using USACE.HEC.Consequences;
 using USACE.HEC.Geography;
+using Geospatial;
+using USACE.HEC.Results;
 
 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 };
+    Geospatial.Utilities.InitializeGDAL();
+
+
+    // 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, Geospatial.Utilities.StructureFieldTypes);
     int count = 0;
-    var watch = System.Diagnostics.Stopwatch.StartNew();
+
     await sp.Process(boundingBox2, (IConsequencesReceptor s) => {
       //Console.WriteLine(((Structure)s).Name);
+      Result res = USACE.HEC.Results.Utilities.ConsequenceReceptorToResult<Structure>(s);
+      c.Write(res);
       count++;
     });
-    watch.Stop();
-    var elapsedMs = watch.ElapsedMilliseconds;
-
     Console.WriteLine(count);
-    Console.WriteLine("Time elapsed: " + elapsedMs.ToString() + " milliseconds");
-    Console.Read(); 
+    Console.Read();
+
+    //Read();
+  }
+
+  public static void Read()
+  {
+
+    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;
+    SpatialProcessor reader = new SpatialProcessor(path);
+    reader.Process<Structure>((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>

ScratchPaper/generated/layer_name.prj

@@ -0,0 +1 @@
+PROJCS["NAD_1983_California_Teale_Albers",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["Albers"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",-4000000.0],PARAMETER["Central_Meridian",-120.0],PARAMETER["Standard_Parallel_1",34.0],PARAMETER["Standard_Parallel_2",40.5],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]