Changeset 4793 for palm/trunk/SCRIPTS/palm_csd_files/palm_csd_canopy_generator.py

Timestamp:

Nov 24, 2020 9:48:21 AM (3 years ago)

Author:

maronga

Message:

allow for overhanging vegetation in palm_csd

File:

• palm/trunk/SCRIPTS/palm_csd_files/palm_csd_canopy_generator.py (modified) (11 diffs)

palm/trunk/SCRIPTS/palm_csd_files/palm_csd_canopy_generator.py

@@ -20 +20 @@
 # Current revisions:
 # -----------------
-# 
-# 
+#
+#
 # Former revisions:
 # -----------------
 # $Id: palm_csd_canopy_generator.py 3773 2019-03-01 08:56:57Z maronga $
+# Bugfix: do not allow BAD < 0, revised treatment of tree trunks
+# 
+# 3773 2019-03-01 08:56:57Z maronga
 # Added support for tree shape
 # 
@@ -46 +49 @@
 import math
 import scipy.integrate as integrate
+import scipy.ndimage as nd
 
 def generate_single_tree_lad(x,y,dz,max_tree_height,max_patch_height,tree_type,tree_height,tree_dia,trunk_dia,tree_lai,season,fill):
@@ -64 +68 @@
    ids = np.ones((len(zlad),len(y),len(x)))
    ids[:,:,:] = fill 
-   
+ 
+   types = np.ones((len(zlad),len(y),len(x)))
+   types[:,:,:] = np.byte(-127)  
 
 #  Calculating the number of trees in the arrays and a boolean array storing the location of trees which is used for convenience in the following loop environment   
@@ -107 +113 @@
                   bad[0:lad_loc_nz,ind_l_y:ind_r_y+1,ind_l_x:ind_r_x+1] = np.where(bad_loc[0:lad_loc_nz,out_l_y:out_r_y+1,out_l_x:out_r_x+1] != fill,bad_loc[0:lad_loc_nz,out_l_y:out_r_y+1,out_l_x:out_r_x+1],bad[0:lad_loc_nz,ind_l_y:ind_r_y+1,ind_l_x:ind_r_x+1]) 
                   ids[0:lad_loc_nz,ind_l_y:ind_r_y+1,ind_l_x:ind_r_x+1] = np.where(lad_loc[0:lad_loc_nz,out_l_y:out_r_y+1,out_l_x:out_r_x+1] != fill,tree_id_counter,ids[0:lad_loc_nz,ind_l_y:ind_r_y+1,ind_l_x:ind_r_x+1]) 
+                  types[0:lad_loc_nz,ind_l_y:ind_r_y+1,ind_l_x:ind_r_x+1] = np.where(lad_loc[0:lad_loc_nz,out_l_y:out_r_y+1,out_l_x:out_r_x+1] != fill,np.byte(tree_type[j,i]),types[0:lad_loc_nz,ind_l_y:ind_r_y+1,ind_l_x:ind_r_x+1])                   
 
 
@@ -116 +123 @@
 
                del lad_loc, x_loc, y_loc, z_loc, status
-   return lad, bad, ids, zlad
+   return lad, bad, ids, types, zlad
 
 
@@ -230 +237 @@
    fillvalues = {
    "tree_data": float(-9999.0),
+   "tree_type": np.byte(-127)
    }
    
    
 #  Check for missing data in the input and set default values if needed
-   if ( tree_type == fillvalues["tree_data"] ):
+   if ( tree_type == fillvalues["tree_type"] ):
       tree_type = int(0)
    else:
@@ -295 +303 @@
    nx = int(tree_dia / dx) + 2
    nz = int(tree_height / dz) + 2
-    
+   
+
+ 
 #  Add one grid point if diameter is an odd value    
    if ( (tree_dia % 2.0) != 0.0 ):
       nx = nx + 1
-
 
 #  Create local domain of the tree's LAD
@@ -459 +468 @@
 
 #  Create BAD array and populate. TODO: revise as low LAD inside the foliage does not result in low BAD values.
-   bad_loc = np.where(lad_loc != fillvalues["tree_data"],(1.0 - lad_loc)*0.1,lad_loc)
+   bad_loc = np.where(lad_loc != fillvalues["tree_data"],(1.0 - (lad_loc/(max(lad_loc.flatten())+0.01)))*0.1,lad_loc)
    
 
@@ -469 +478 @@
             if ( z[k] <= crown_center ):
                r_test = np.sqrt( (x[i] - tree_location_x)**(2) + (y[j] - tree_location_y)**(2)  )
-               if ( r_test <= radius ):
-                  bad_loc[k,j,i] = 1.0
-               if ( (r_test == 0.0) & (trunk_dia <= dx) ):
-                   bad_loc[k,j,i] = radius**(2) * 3.14159265359
+               if ( r_test == 0.0 ):
+                  if ( trunk_dia <= dx ):
+                     bad_loc[k,j,i] = radius**2 * 3.14159265359
+                  else:
+                     # WORKAROUND: divide remaining circle area over the 8 surrounding valid_pixels
+                     bad_loc[k,j-1:j+2,i-1:i+2] = radius**2 * 3.14159265359 / 8.0
+                     # for the central pixel fill the pixel
+                     bad_loc[k,j,i] = dx**2
+               #elif ( r_test <= radius ):
+                  # TODO: calculate circle segment of grid points cut by the grid
+ 
 
    return lad_loc, bad_loc, x, y, z, 0
 
 
-def process_patch(dz,patch_height,max_height_lad,patch_lai,alpha,beta):
+def process_patch(dz,patch_height,vegetation_type,max_height_lad,patch_lai,alpha,beta):
    
 #  Define fill values
    fillvalues = {
    "tree_data": float(-9999.0),
-   "pch_index": int(-9999),
+   "tree_id":   int(-9999),    
+   "tree_type": np.byte(-127),  
+   "pch_index": int(-9999)
    }
 
+   
    phdz = patch_height[:,:] / dz
    pch_index = np.where( (patch_height[:,:] != fillvalues["tree_data"]),phdz.astype(int)+1,int(-1))    
@@ -520 +539 @@
                 lad_loc[k,j,i] = 0.5 * ( pre_lad[k-1] + pre_lad[k] )
 
-   return lad_loc, nz, 0
+
+   patch_id_2d = np.where(lad_loc[0,:,:] != fillvalues["tree_data"], 1, 0)
+   patch_id = np.where(lad_loc != fillvalues["tree_data"], 1, 0)
+   patch_type = np.where(lad_loc != fillvalues["tree_data"], np.byte(1), np.byte(-127))
+
+   num_ids = nd.label(patch_id_2d, output=patch_id_2d)
+ 
+   for k in range(0,nz):
+   
+      vegetation_type_int = vegetation_type.astype(int)
+      patch_id[k,:,:]   = np.where((patch_id_2d != 0) & (lad_loc[k,:,:] != fillvalues["tree_data"]), patch_id_2d, fillvalues["tree_id"])
+      patch_type[k,:,:] = np.where((patch_id_2d != 0) & (lad_loc[k,:,:] != fillvalues["tree_data"]), np.byte(100+vegetation_type), fillvalues["tree_type"])
+
+   return lad_loc, patch_id, patch_type, nz, 0
 
 
@@ -559 +591 @@
        self.bad_scale = bad_scale
        self.dbh = dbh    
+