Changeset 4793 for palm/trunk/SCRIPTS/palm_csd

Timestamp:

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

Author:

maronga

Message:

allow for overhanging vegetation in palm_csd

File:

• palm/trunk/SCRIPTS/palm_csd (modified) (15 diffs)

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