gradient search fix and nc_masked_kml

.github/workflows/increment-version.yml

@@ -2,7 +2,7 @@ name: Increment Version on Dev
 
 on:
   push:
-    branches: [ dev ]
+    # branches: [ dev ] we can add back later if needed
   workflow_dispatch:
 
 jobs:

docs/source/about/statement_of_need.rst

@@ -26,8 +26,4 @@ ForeFire is intended for:
 
 *   **Wildfire Researchers:** Investigating fire behavior, model development, sensitivity analyses, fire-atmosphere interactions.
 *   **Operational Agencies:** Use in forecasting systems (potentially via ensemble simulations), risk assessment, and post-fire analysis.
-*   **Students and Educators:** A tool for learning about wildfire dynamics and simulation techniques.
-
-**Context:**
-
-ForeFire aims to complement and extend the capabilities of other existing wildfire simulation tools by focusing on high-performance computing, advanced physics coupling, and providing a flexible open-source C++ core. *(Mention specific limitations of other tools that ForeFire addresses, if appropriate, and potentially cite 1-2 key alternative tools for comparison).*
+*   **Students and Educators:** A tool for learning about wildfire dynamics and simulation techniques.

src/GradientDataLayer.h

@@ -103,11 +103,12 @@ T GradientDataLayer<T>::getValueAt(FFPoint loc, const double& time) {
     // Compute the gradient in each direction and find the maximum.
     T maxGradient = 0.0;
     for (size_t i = 0; i < directions.size(); ++i) {
-        FFPoint nextLoc = loc + (10 * directions[i]);
+		double dxGrad = getDx();
+        FFPoint nextLoc = loc + (dxGrad * directions[i]);
         T neighborValue = parent->getValueAt(nextLoc, time);
         // Since the displacement magnitude is dx in every case (after normalization),
         // the gradient is computed as the difference divided by dx.
-        T gradient = (neighborValue - currentValue) / 10;
+        T gradient = (neighborValue - currentValue) / dxGrad;
         if (gradient > maxGradient) {
             maxGradient = gradient;
         }

src/include/Version.h

@@ -1,3 +1,3 @@
 #pragma once
-const char* ff_version = "v2.1.109";
+const char* ff_version = "v2.1.114";
 

tools/preprocessing/extract_from_eu_land_cover.py

@@ -28,6 +28,9 @@
 
 from rasterio.features import rasterize
 from shapely.geometry import mapping
+
+import xarray as xr
+
 
 attribute_widths_road_edge_full = {
     'secondary': 0.8,
@@ -552,7 +555,7 @@ def fake_fuel(legend_file_path, WSEN, LBRT,output_dir,fuel_resolution = 10):
 
 
 
-def rasterize_kml(kml_path, ref_tif, output_path, default_value=1, imbounds=None):
+def rasterize_kml(kml_path, ref_tif, output_path, default_value=62, imbounds=None):
     print(f"Rasterizing KML {kml_path} to {output_path}")
 
     # Load reference GeoTIFF for spatial reference
@@ -590,6 +593,57 @@ def rasterize_kml(kml_path, ref_tif, output_path, default_value=1, imbounds=None
                        crs=ref_crs, transform=transform) as dst:
         dst.write(raster, 1)
 
+def nc_masked_kml(ref_nc, path_kml, field='fuel', default_value=1, imbounds=None):
+    import fiona
+    fiona.drvsupport.supported_drivers['KML'] = 'rw'
+
+    if field == 'fuel':
+        ny = ref_nc[field].sizes['fy']
+        nx = ref_nc[field].sizes['fx']
+
+    else:
+        ny = ref_nc.sizes['ny']
+        nx = ref_nc.sizes['nx']
+
+    attrs = ref_nc['domain'].attrs
+
+    west, south, east, north = map(float, attrs['BBoxWSEN'].split(','))
+    if imbounds is not None:
+        west, south, east, north = imbounds
+
+    resx = (east - west) / nx
+    resy = (north - south) / ny
+
+    transform = from_origin(west, north, resx, resy)
+
+
+    gdf = gpd.read_file( path_kml, driver='KML')
+
+
+    shapes = []
+    for _, feature in gdf.iterrows():
+        name = feature.get('Name', '')
+        try:
+            val = int(name) if name.isdigit() else default_value
+        except Exception:
+            val = default_value
+        geom = feature.geometry
+        shapes.append((geom, val))
+
+    rasterized = rasterize(
+            shapes,
+            out_shape=(ny, nx),
+            transform=transform,
+            fill=0,
+            dtype='int16')
+
+    fieldnc = ref_nc[field][0,0,::-1,:] #for fuels it is upside down
+    fieldnc_updated = xr.where(rasterized == 1, 62, fieldnc)
+
+    ref_nc[field][0,0,:,:] = fieldnc_updated[::-1,:]
+
+    return ref_nc
+
 def landcover_roads_to_fuel(S2GLC_tif,legend_file_path, WSEN, LBRT,output_dir,fuel_modifier=None, fuel_resolution = 10, no_fuel_code = 62):
 
     fuel_indices_origin = f"{output_dir}/fuel_indices_S2GLC.tif"