From 67017d092329a5bfaa8a5b88a8d391ca69d5ad11 Mon Sep 17 00:00:00 2001 From: Tobey Carman Date: Wed, 21 Jul 2021 13:16:20 -0800 Subject: [PATCH 1/2] Change to using km-2 for "Area of Burn" after... ...reviewing code and discussing with H. Genet. --- scripts/create_region_input.py | 10 +++++----- src/WildFire.cpp | 1 + 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/scripts/create_region_input.py b/scripts/create_region_input.py index 57c721957..3ca0bf0aa 100755 --- a/scripts/create_region_input.py +++ b/scripts/create_region_input.py @@ -1561,9 +1561,9 @@ def fill_explicit_fire_file(startyr, yrs, xo, yo, xs, ys, out_dir, of_name, tiff # Open the temporary file and read the transformation matrix. Use this # matrix with the rasterize function to geo-reference the features. rmeta = rasterio.open(tmpFile) - # NOTE THAT THE AREA HERE IS IN m^2, not km^2 !!! - shapes = [(geom, value) for geom, value in zip(this_years_fires.geometry, this_years_fires.Shape_Area)] - + # NOTE: Area of Burn stored in the shape file from the BLM is + # expressed in m^2. We convert it here to km^2 as required by dvmdostem. + shapes = [(geom, value) for geom, value in zip(this_years_fires.geometry, this_years_fires.Shape_Area/1000.0/1000.0)] aob = features.rasterize(shapes, out_shape=(ys,xs), transform=rmeta.transform) severity = np.greater(aob,0) * 2 # Sets any pixel with area of burn > 0 to severity of 2 jday = np.greater(aob,0) * 212 @@ -1625,7 +1625,8 @@ def fill_explicit_fire_file(startyr, yrs, xo, yo, xs, ys, out_dir, of_name, tiff print("Building scarnum2area_map based on full spatial domain...") # make a set of all the FIREIDs in the file, then map to the count of - # the pixels (1km pixels) that share FIREID + # the pixels (1km pixels) that share FIREID. The result is that the + # area is expressed here as km^2. scarnum2area_map = {} with netCDF4.Dataset(tmp_fs_spatial_full) as ds: fs_data = ds.variables['Band2'][:] @@ -1676,7 +1677,6 @@ def lookup_AOB(x): # correct. I.e. 2 pixels (2 km-2) of a 100 pixel (100 km-2) burn # fall in the users selected area --> AOB will be 100 for each pixel. aob = lookup_AOB_v(fs_d) - aob = aob * 1000 * 1000 # Convert from square kilometers to square meters # Now get the burn severity with netCDF4.Dataset(tmp_bs_file, 'r') as bs_ds: diff --git a/src/WildFire.cpp b/src/WildFire.cpp index b8ce89a5e..8c6df35b0 100644 --- a/src/WildFire.cpp +++ b/src/WildFire.cpp @@ -619,6 +619,7 @@ double WildFire::getBurnOrgSoilthick(const int year) { folb = 0.01*((21.5-6.1)/21.5); } else { // boreal forest: Genet et al.2013; if(this->slope<2){ + // AOB in km-2, see coefficient and paper. folb = 0.1276966713-0.0319397467*this->slope+0.0020914862*this->exp_jday_of_burn[year]+0.0127016155* log(this->exp_area_of_burn[year]); } else { folb = -0.2758306315+0.0117609336*this->slope-0.0744057680*cos(this->asp * 3.14159265 / 180 ) +0.0260221684*edall->m_soid.tshlw+0.0011413114*this->exp_jday_of_burn[year]+0.0336302905*log (this->exp_area_of_burn[year]); From bb6a162165216786b342ceee003c116840ceacdd Mon Sep 17 00:00:00 2001 From: Tobey Carman Date: Wed, 21 Jul 2021 13:19:16 -0800 Subject: [PATCH 2/2] Fix bug with severity, change default severity to '3'... ... for historic observed fires After discussing with H. Genet and reviewing modeled fire severity. Bug related to reading/writing masked arrays from netCDF files. Without setting missing_value attribute on the variable, the variable ends up being written data from the mask, (i.e. all ones) and not the underlying data. --- scripts/create_region_input.py | 3 ++- src/WildFire.cpp | 1 + 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/scripts/create_region_input.py b/scripts/create_region_input.py index 3ca0bf0aa..5b13624c6 100755 --- a/scripts/create_region_input.py +++ b/scripts/create_region_input.py @@ -1565,7 +1565,7 @@ def fill_explicit_fire_file(startyr, yrs, xo, yo, xs, ys, out_dir, of_name, tiff # expressed in m^2. We convert it here to km^2 as required by dvmdostem. shapes = [(geom, value) for geom, value in zip(this_years_fires.geometry, this_years_fires.Shape_Area/1000.0/1000.0)] aob = features.rasterize(shapes, out_shape=(ys,xs), transform=rmeta.transform) - severity = np.greater(aob,0) * 2 # Sets any pixel with area of burn > 0 to severity of 2 + severity = np.greater(aob,0) * 3 # default to burn severity of 3 jday = np.greater(aob,0) * 212 mask = np.greater(aob,0) @@ -1685,6 +1685,7 @@ def lookup_AOB(x): # Write AOB and Burn Severity to output file... with netCDF4.Dataset(of_name, 'a') as ds: ds.variables['exp_area_of_burn'][iy,:] = aob + ds.variables['exp_fire_severity'].missing_value = 255 ds.variables['exp_fire_severity'][iy,:] = severity ds.variables['exp_jday_of_burn'][iy,:] = np.greater(aob,0) * 212 # July 31 2021 ds.variables['exp_burn_mask'][iy,:] = np.logical_not(aob.mask) diff --git a/src/WildFire.cpp b/src/WildFire.cpp index 8c6df35b0..74e252ade 100644 --- a/src/WildFire.cpp +++ b/src/WildFire.cpp @@ -578,6 +578,7 @@ double WildFire::getBurnOrgSoilthick(const int year) { double burn_thickness = 0.0; // For now, severity class is based on ALFRESCO: + // NOTE, this may be incorrect, it looks like ALFRESCO might use 1 for no burn. // 0 - no burning // 1 - low // 2 - moderate