Changeset 4793 for palm/trunk/SCRIPTS

Timestamp:

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

Author:

maronga

Message:

allow for overhanging vegetation in palm_csd

Location:

palm/trunk/SCRIPTS

Files:

• .csd.config.default (modified) (2 diffs)
• palm_csd (modified) (15 diffs)
• palm_csd_files/palm_csd_canopy_generator.py (modified) (11 diffs)

palm/trunk/SCRIPTS/.csd.config.default

@@ -212 +212 @@
 vegetation_on_roofs = False
 street_trees = True
+overhanging_trees = False
 
 [domain_N02]
@@ -228 +229 @@
 vegetation_on_roofs = False
 street_trees = True
+overhanging_trees = False

palm/trunk/SCRIPTS/palm_csd

@@ -25 +25 @@
 # -----------------
 # $Id$
+# Added output of tree_type (species or classification). Added tree_id for
+# vegetation patches. Allow for overlapping vegetation with buildings
+#
+# 4749 2020-10-16 11:21:27Z maronga
 # Bugfix: processing of patch height faulty for multiple domains
 #
@@ -164 +168 @@
    global domain_street_trees
    global domain_canopy_patches
+   global domain_overhanging_trees
 
 #  Definition of input data parameters
@@ -250 +255 @@
    domain_street_trees = []
    domain_canopy_patches = []
+   domain_overhanging_trees = []
 
    zt_min = 0.0
@@ -349 +355 @@
          domain_green_roofs.append(config.getboolean(read_tmp, 'vegetation_on_roofs')) 
          domain_street_trees.append(config.getboolean(read_tmp, 'street_trees')) 
-         
+         domain_overhanging_trees.append(config.getboolean(read_tmp, 'overhanging_trees')) 
+                 
       if ( read_tmp.split("_")[0] == 'input' ):
          input_names.append(read_tmp.split("_")[1])
@@ -417 +424 @@
    "vegetation_pars": "f4",
    "tree_data": "f4",
-   "tree_type": "b",
+   "tree_id": "i", 
+   "tree_type": "f4",
+   "tree_types": "b",  
    "nbuilding_pars": "i",
    "nvegetation_pars": "i",
@@ -447 +456 @@
    "vegetation_pars": float(-9999.0),
    "tree_data": float(-9999.0),
-   "tree_type": np.byte(-127)
+   "tree_id": int(-9999),
+   "tree_type": float(-9999.0),
+   "tree_types": np.byte(-127)
    }
 
@@ -473 +484 @@
    "buildings_pars": float(-9999.0),
    "tree_data": float(-9999.0),
-   "tree_type": 0,
+   "tree_type": np.byte(-127),
    "vegetation_pars": float(-9999.0)   
    }
@@ -963 +974 @@
       vegetation_pars =  nc_read_from_file_2d_all(filename[i], 'vegetation_pars')
 
-      lai = np.copy(vegetation_pars[1,:,:])
+      lai = nc_read_from_file_2d(input_file_lai[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])   
 
       x = nc_read_from_file_2d_all(filename[i], 'x')
@@ -991 +1002 @@
       tree_trunk_diameter = np.where( (tree_trunk_diameter == 0.0) | (tree_trunk_diameter == -1.0) ,fillvalues["tree_data"],tree_trunk_diameter)
       tree_type = np.where( (tree_type == 0.0) | (tree_type == -1.0) ,fillvalues["tree_data"],tree_type)
-      patch_height = np.where( patch_height == -1.0 ,fillvalues["tree_data"],patch_height)     
+      
+#     For vegetation pixel with missing height information (-1), set default patch height   
+      patch_height = np.where( patch_height == -1.0 ,settings_patch_height_default,patch_height)     
 
 #     Convert trunk diameter from cm to m
@@ -1005 +1018 @@
 
 #     For zero-height patches, set vegetation_type to short grass and remove these pixels from the patch height array
-      vegetation_type = np.where( (patch_height == 0.0) & ( (vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18) ),3,vegetation_type)
-      patch_type = np.where( (patch_height == 0.0) & ( (vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18) ),fillvalues["tree_data"],patch_height)   
-      
+      vegetation_type = np.where( (patch_height == 0.0) & ( (vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18) ),3,vegetation_type) 
+      patch_height = np.where( (patch_height == 0.0) & ( (vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18) ),fillvalues["tree_data"],patch_height)       
 
       max_tree_height = max(tree_height.flatten())
@@ -1017 +1029 @@
       if ( (max_tree_height != fillvalues["tree_data"]) | (max_patch_height != fillvalues["tree_data"]) ):
 
-         lad, bad, tree_ids, zlad = generate_single_tree_lad(x,y,domain_dz[i],max_tree_height,max_patch_height,tree_type,tree_height,tree_crown_diameter,tree_trunk_diameter,lai,settings_season,fillvalues["tree_data"])
-  
+         lad, bad, tree_ids, tree_types, zlad = generate_single_tree_lad(x,y,domain_dz[i],max_tree_height,max_patch_height,tree_type,tree_height,tree_crown_diameter,tree_trunk_diameter,lai,settings_season,fillvalues["tree_data"])
   
 #        Remove LAD volumes that are inside buildings
-         buildings_2d = nc_read_from_file_2d_all(filename[i], 'buildings_2d') 
-         for k in range(0,len(zlad)-1):
+         if not domain_overhanging_trees[i]:
+            buildings_2d = nc_read_from_file_2d_all(filename[i], 'buildings_2d') 
+            for k in range(0,len(zlad)-1):
   
-            lad[k,:,:] = np.where(buildings_2d == fillvalues["buildings_2d"],lad[k,:,:],fillvalues["tree_data"])
-            bad[k,:,:] = np.where(buildings_2d == fillvalues["buildings_2d"],bad[k,:,:],fillvalues["tree_data"])
-            tree_ids[k,:,:] = np.where(buildings_2d == fillvalues["buildings_2d"],tree_ids[k,:,:],fillvalues["tree_data"])      
-
-         del buildings_2d
+               lad[k,:,:] = np.where(buildings_2d == fillvalues["buildings_2d"],lad[k,:,:],fillvalues["tree_data"])
+               bad[k,:,:] = np.where(buildings_2d == fillvalues["buildings_2d"],bad[k,:,:],fillvalues["tree_data"])
+               tree_ids[k,:,:] = np.where(buildings_2d == fillvalues["buildings_2d"],tree_ids[k,:,:],fillvalues["tree_id"])      
+
+            del buildings_2d
   
+ 
          nc_write_dimension(filename[i], 'zlad', zlad, datatypes["tree_data"]) 
+
          nc_write_to_file_3d(filename[i], 'lad', lad, datatypes["tree_data"],'zlad','y','x',fillvalues["tree_data"]) 
          nc_write_attribute(filename[i], 'lad', 'long_name', 'leaf area density')
@@ -1037 +1051 @@
          nc_write_attribute(filename[i], 'lad', 'coordinates', 'E_UTM N_UTM lon lat') 
          nc_write_attribute(filename[i], 'lad', 'grid_mapping', 'E_UTM N_UTM lon lat') 
-  
+ 
          nc_write_to_file_3d(filename[i], 'bad', bad, datatypes["tree_data"],'zlad','y','x',fillvalues["tree_data"]) 
          nc_write_attribute(filename[i], 'bad', 'long_name', 'basal area density')
@@ -1045 +1059 @@
          nc_write_attribute(filename[i], 'bad', 'grid_mapping', 'E_UTM N_UTM lon lat')
             
-         nc_write_to_file_3d(filename[i], 'tree_id', tree_ids, datatypes["tree_data"],'zlad','y','x',fillvalues["tree_data"]) 
+         nc_write_to_file_3d(filename[i], 'tree_id', tree_ids, datatypes["tree_id"],'zlad','y','x',fillvalues["tree_id"]) 
          nc_write_attribute(filename[i], 'tree_id', 'long_name', 'tree id')
          nc_write_attribute(filename[i], 'tree_id', 'units', '')
@@ -1051 +1065 @@
          nc_write_attribute(filename[i], 'tree_id', 'coordinates', 'E_UTM N_UTM lon lat') 
          nc_write_attribute(filename[i], 'tree_id', 'grid_mapping', 'E_UTM N_UTM lon lat')           
+
+         nc_write_to_file_3d(filename[i], 'tree_type', tree_types, datatypes["tree_types"],'zlad','y','x',fillvalues["tree_types"]) 
+         nc_write_attribute(filename[i], 'tree_type', 'long_name', 'tree type')
+         nc_write_attribute(filename[i], 'tree_type', 'units', '')
+         nc_write_attribute(filename[i], 'tree_type', 'res_orig', domain_px[i])  
+         nc_write_attribute(filename[i], 'tree_type', 'coordinates', 'E_UTM N_UTM lon lat') 
+         nc_write_attribute(filename[i], 'tree_type', 'grid_mapping', 'E_UTM N_UTM lon lat')   
             
-         del lai, lad, bad, tree_ids, zlad
-
+         del lai, lad, bad, tree_ids, tree_types, zlad
+
+      else:
+         print('No street trees generated in domain ' + str(i))
+         
       del vegetation_pars, tree_height, tree_crown_diameter, tree_trunk_diameter, tree_type, patch_height, x, y
 
@@ -1061 +1085 @@
    if domain_canopy_patches[i]: 
 
-#     Load vegetation_type and lad array (at level z = 0) for re-processing
-      vegetation_type = nc_read_from_file_2d_all(filename[i], 'vegetation_type') 
-      lad = nc_read_from_file_3d_all(filename[i], 'lad') 
-      zlad = nc_read_from_file_1d_all(filename[i], 'zlad') 
       patch_height = nc_read_from_file_2d(input_file_patch_height[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) 
-      vegetation_pars = nc_read_from_file_3d_all(filename[i], 'vegetation_pars') 
-      lai = vegetation_pars[1,:,:]
-     
-
-#     Determine all pixels that do not already have an LAD but which are high vegetation to a dummy value of 1.0 and remove all other pixels
-      lai_high = np.where( (lad[0,:,:] == fillvalues["tree_data"]) & ( ( (vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18) ) & ( (patch_height == fillvalues["tree_data"]) | (patch_height >= domain_dz[i])) ),1.0,fillvalues["tree_data"])
-
-
-#     Treat all pixels where short grass is defined, but where a patch_height >= dz is found, as high vegetation (often the case in backyards)
-      lai_high = np.where( (lai_high == fillvalues["tree_data"]) & (patch_height >= domain_dz[i]) & (vegetation_type == 3), 1.0, lai_high)
+
+#     For vegetation pixel with missing height information (-1), set default patch height      
+      patch_height = np.where( patch_height == -1.0 ,settings_patch_height_default,patch_height)
+      
+      max_patch_height = max(patch_height.flatten())
+      if ( max_patch_height != fillvalues["tree_data"] ):
+#        Call canopy generator for single trees only if there is any tree height available in the domain.
+#        This does not guarantee that there are street trees that can be processed. This is checked in the
+#        canopy generator.
+
+#        Load vegetation_type and lad array (at level z = 0) for re-processing
+         vegetation_type = nc_read_from_file_2d_all(filename[i], 'vegetation_type') 
+         lad = nc_read_from_file_3d_all(filename[i], 'lad') 
+         tree_id = nc_read_from_file_3d_all(filename[i], 'tree_id') 
+         tree_types = nc_read_from_file_3d_all(filename[i], 'tree_type') 
+         zlad = nc_read_from_file_1d_all(filename[i], 'zlad') 
+         vegetation_pars = nc_read_from_file_3d_all(filename[i], 'vegetation_pars') 
+         lai = nc_read_from_file_2d(input_file_lai[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])   
+    
+         nc_write_to_file_2d(filename[i], 'patch_height', patch_height, datatypes["tree_data"],'y','x',fillvalues["tree_data"]) 
+         nc_write_attribute(filename[i], 'patch_height', 'long_name', 'tree type')
+         nc_write_attribute(filename[i], 'patch_height', 'units', '')
+         nc_write_attribute(filename[i], 'patch_height', 'res_orig', domain_px[i])  
+         nc_write_attribute(filename[i], 'patch_height', 'coordinates', 'E_UTM N_UTM lon lat') 
+         nc_write_attribute(filename[i], 'patch_height', 'grid_mapping', 'E_UTM N_UTM lon lat')      
+
+#        Determine all pixels that do not already have an LAD but which are high vegetation to a dummy value of 1.0 and remove all other pixels
+         lai_high = np.where( (lad[0,:,:] == fillvalues["tree_data"]) & ( ( (vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18) ) & ( (patch_height == fillvalues["tree_data"]) | (patch_height >= domain_dz[i])) ),1.0,fillvalues["tree_data"])
+
+
+#        Treat all pixels where short grass is defined, but where a patch_height >= dz is found, as high vegetation (often the case in backyards)
+         lai_high = np.where( (lai_high == fillvalues["tree_data"]) & (patch_height >= domain_dz[i]) & (vegetation_type == 3), 1.0, lai_high)
                           
-#     Now, assign either the default LAI for high vegetation or keep 1.0 from the lai_high array.
-      lai_high = np.where( (lai_high != fillvalues["tree_data"]) & (lai == fillvalues["tree_data"]), settings_lai_high_default, lai_high)
-
-#     If LAI values are available in the LAI array, write them on the lai_high array
-      lai_high = np.where( (lai_high != fillvalues["tree_data"]) & (lai != fillvalues["tree_data"]), lai, lai_high)
-
-
-#     Define a patch height wherever it is missing, but where a high vegetation LAI was set
-      patch_height = np.where ( (lai_high != fillvalues["tree_data"]) & (patch_height == fillvalues["tree_data"]), settings_patch_height_default, patch_height)
-
-
-#     Remove pixels where street trees were already set
-      patch_height = np.where ( (lad[0,:,:] != fillvalues["tree_data"]), fillvalues["tree_data"], patch_height)
-
-#     Remove patch heights that have no lai_high value
-      patch_height = np.where ( (lai_high == fillvalues["tree_data"]), fillvalues["tree_data"], patch_height)
-
-#     For missing LAI values, set either the high vegetation default or the low vegetation default    
-      lai_high = np.where( (patch_height > 2.0) & (patch_height != fillvalues["tree_data"]) & (lai_high == fillvalues["tree_data"]),settings_lai_high_default,lai_high)
-      lai_high = np.where( (patch_height <= 2.0) & (patch_height != fillvalues["tree_data"]) & (lai_high == fillvalues["tree_data"]),settings_lai_low_default,lai_high)
+#        If overhanging trees are allowed, assign pixels with patch_height > dz that are not included in vegetation_type to high
+         if domain_overhanging_trees[i]:
+            lai_high = np.where( (lai_high == fillvalues["tree_data"]) & (patch_height >= domain_dz[i]), 1.0, lai_high)
+
+#        Now, assign either the default LAI for high vegetation or keep 1.0 from the lai_high array.
+         lai_high = np.where( (lai_high != fillvalues["tree_data"]) & (lai == fillvalues["tree_data"]), settings_lai_high_default, lai_high)
+
+#        If LAI values are available in the LAI array, write them on the lai_high array
+         lai_high = np.where( (lai_high != fillvalues["tree_data"]) & (lai != fillvalues["tree_data"]), lai, lai_high)
+
+
+#        Define a patch height wherever it is missing, but where a high vegetation LAI was set
+         patch_height = np.where ( (lai_high != fillvalues["tree_data"]) & (patch_height == fillvalues["tree_data"]), settings_patch_height_default, patch_height)
+
+
+#        Remove pixels where street trees were already se
+         patch_height = np.where ( (lad[0,:,:] != fillvalues["tree_data"]), fillvalues["tree_data"], patch_height)
+
+#        Remove patch heights that have no lai_high value
+         patch_height = np.where ( (lai_high == fillvalues["tree_data"]), fillvalues["tree_data"], patch_height)
+
+
+
+
+#        For missing LAI values, set either the high vegetation default or the low vegetation default    
+         lai_high = np.where( (patch_height > 2.0) & (patch_height != fillvalues["tree_data"]) & (lai_high == fillvalues["tree_data"]),settings_lai_high_default,lai_high)
+         lai_high = np.where( (patch_height <= 2.0) & (patch_height != fillvalues["tree_data"]) & (lai_high == fillvalues["tree_data"]),settings_lai_low_default,lai_high)
          
-      if ( max(patch_height.flatten()) >= (2.0 * domain_dz[i]) ):
-         print("    start calculating LAD (this might take some time)")       
+         if ( max(patch_height.flatten()) >= (2.0 * domain_dz[i]) ):
+            print("    start calculating LAD (this might take some time)")       
 
          
-         lad_patch, patch_nz, status = process_patch(domain_dz[i],patch_height,max(zlad),lai_high,settings_lai_alpha,settings_lai_beta)
-
-#        2D loop in order to avoid memory problems with large arrays
-         for iii in range(0,domain_nx[i]):
-            for jj in range(0,domain_ny[i]):            
-
-               lad[0:patch_nz+1,jj,iii] = np.where( (lad[0:patch_nz+1,jj,iii] == fillvalues["tree_data"]),lad_patch[0:patch_nz+1,jj,iii], lad[0:patch_nz+1,jj,iii] )
-
-#     Remove high vegetation wherever it is replaced by a leaf area density. This should effectively remove all high vegetation pixels
-      vegetation_type = np.where((lad[0,:,:] != fillvalues["tree_data"]) & (vegetation_type != fillvalues["vegetation_type"]),settings_veg_type_below_trees,vegetation_type)
+            lad_patch, patch_id, patch_types, patch_nz, status = process_patch(domain_dz[i],patch_height,vegetation_type,max(zlad),lai_high,settings_lai_alpha,settings_lai_beta)
+
+#           Increase the patch ids by the maximum id for street trees to avoid double ids
+            max_id = max(tree_id.flatten())
+            patch_id = np.where( patch_id != fillvalues["tree_id"], patch_id+max_id, patch_id)
+
+#           2D loop in order to avoid memory problems with large arrays
+            for iii in range(0,domain_nx[i]):
+               for jj in range(0,domain_ny[i]):            
+
+                  tree_id[0:patch_nz+1,jj,iii] = np.where( (lad[0:patch_nz+1,jj,iii] == fillvalues["tree_data"]),patch_id[0:patch_nz+1,jj,iii], tree_id[0:patch_nz+1,jj,iii] )
+                  tree_types[0:patch_nz+1,jj,iii] = np.where( (lad[0:patch_nz+1,jj,iii] == fillvalues["tree_data"]),patch_types[0:patch_nz+1,jj,iii], tree_types[0:patch_nz+1,jj,iii] )                  
+                  lad[0:patch_nz+1,jj,iii]     = np.where( (lad[0:patch_nz+1,jj,iii] == fillvalues["tree_data"]),lad_patch[0:patch_nz+1,jj,iii], lad[0:patch_nz+1,jj,iii] )
+
+#        Remove high vegetation wherever it is replaced by a leaf area density. This should effectively remove all high vegetation pixels
+         vegetation_type = np.where((lad[0,:,:] != fillvalues["tree_data"]) & (vegetation_type != fillvalues["vegetation_type"]),settings_veg_type_below_trees,vegetation_type)
               
-#     If desired, remove all high vegetation. TODO: check if this is still necessary
-      if not domain_high_vegetation[i]:
-         vegetation_type = np.where((vegetation_type != fillvalues["vegetation_type"]) & ( (vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18) ),3,vegetation_type)   
+#        If desired, remove all high vegetation. TODO: check if this is still necessary
+         if not domain_high_vegetation[i]:
+            vegetation_type = np.where((vegetation_type != fillvalues["vegetation_type"]) & ( (vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18) ),3,vegetation_type)   
 
   
-#     Set default low LAI for pixels with an LAD (short grass below trees)
-      lai_low = np.where( (lad[0,:,:] == fillvalues["tree_data"]), lai, settings_lai_low_default)
+#        Set default low LAI for pixels with an LAD (short grass below trees)
+         lai_low = np.where( (lad[0,:,:] == fillvalues["tree_data"]), lai, settings_lai_low_default)
     
-#     Fill low vegetation pixels without LAI set or with LAI = 0 with default value
-      lai_low = np.where( ( (lai_low == fillvalues["tree_data"]) | (lai_low == 0.0) ) & (vegetation_type != fillvalues["vegetation_type"] ), settings_lai_low_default, lai_low)
-
-#     Remove lai for pixels that have no vegetation_type
-      lai_low = np.where( ( (vegetation_type != fillvalues["vegetation_type"]) & (vegetation_type != 1) ), lai_low, fillvalues["tree_data"])
+#        Fill low vegetation pixels without LAI set or with LAI = 0 with default value
+         lai_low = np.where( ( (lai_low == fillvalues["tree_data"]) | (lai_low == 0.0) ) & (vegetation_type != fillvalues["vegetation_type"] ), settings_lai_low_default, lai_low)
+
+#        Remove lai for pixels that have no vegetation_type
+         lai_low = np.where( ( (vegetation_type != fillvalues["vegetation_type"]) & (vegetation_type != 1) ), lai_low, fillvalues["tree_data"])
         
-#     Overwrite lai in vegetation_parameters  
-      vegetation_pars[1,:,:] = np.copy(lai_low)
-      nc_overwrite_to_file_3d(filename[i], 'vegetation_pars', vegetation_pars) 
-
-#     Overwrite lad array
-      nc_overwrite_to_file_3d(filename[i], 'lad', lad) 
-           
-      nc_overwrite_to_file_2d(filename[i], 'vegetation_type', vegetation_type)  
-
-      
-      del vegetation_type, lad, lai, patch_height, vegetation_pars, zlad
-
+#        Overwrite lai in vegetation_parameters  
+         vegetation_pars[1,:,:] = np.copy(lai_low)
+         nc_overwrite_to_file_3d(filename[i], 'vegetation_pars', vegetation_pars) 
+
+#        Overwrite lad and id arrays
+         nc_overwrite_to_file_3d(filename[i], 'lad', lad) 
+         nc_overwrite_to_file_3d(filename[i], 'tree_id', tree_id) 
+         nc_overwrite_to_file_3d(filename[i], 'tree_type', tree_types) 
+
+         nc_overwrite_to_file_2d(filename[i], 'vegetation_type', vegetation_type)  
+
+
+         del vegetation_type, lad, lai, patch_height, vegetation_pars, zlad
+      else:
+         print('No tree patches found in domain ' + str(i))
 
 # Final adjustment of vegetation parameters: remove LAI where a bare soil was set

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    
+