Ignore:
Timestamp:
Nov 24, 2020 9:48:21 AM (23 months ago)
Author:
maronga
Message:

allow for overhanging vegetation in palm_csd

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SCRIPTS/palm_csd_files/palm_csd_canopy_generator.py

    r4481 r4793  
    2020# Current revisions:
    2121# -----------------
    22 # 
    23 # 
     22#
     23#
    2424# Former revisions:
    2525# -----------------
    2626# $Id: palm_csd_canopy_generator.py 3773 2019-03-01 08:56:57Z maronga $
     27# Bugfix: do not allow BAD < 0, revised treatment of tree trunks
     28#
     29# 3773 2019-03-01 08:56:57Z maronga
    2730# Added support for tree shape
    2831#
     
    4649import math
    4750import scipy.integrate as integrate
     51import scipy.ndimage as nd
    4852
    4953def 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):
     
    6468   ids = np.ones((len(zlad),len(y),len(x)))
    6569   ids[:,:,:] = fill
    66    
     70 
     71   types = np.ones((len(zlad),len(y),len(x)))
     72   types[:,:,:] = np.byte(-127) 
    6773
    6874#  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   
     
    107113                  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])
    108114                  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])
     115                  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])                   
    109116
    110117
     
    116123
    117124               del lad_loc, x_loc, y_loc, z_loc, status
    118    return lad, bad, ids, zlad
     125   return lad, bad, ids, types, zlad
    119126
    120127
     
    230237   fillvalues = {
    231238   "tree_data": float(-9999.0),
     239   "tree_type": np.byte(-127)
    232240   }
    233241   
    234242   
    235243#  Check for missing data in the input and set default values if needed
    236    if ( tree_type == fillvalues["tree_data"] ):
     244   if ( tree_type == fillvalues["tree_type"] ):
    237245      tree_type = int(0)
    238246   else:
     
    295303   nx = int(tree_dia / dx) + 2
    296304   nz = int(tree_height / dz) + 2
    297    
     305   
     306
     307 
    298308#  Add one grid point if diameter is an odd value   
    299309   if ( (tree_dia % 2.0) != 0.0 ):
    300310      nx = nx + 1
    301 
    302311
    303312#  Create local domain of the tree's LAD
     
    459468
    460469#  Create BAD array and populate. TODO: revise as low LAD inside the foliage does not result in low BAD values.
    461    bad_loc = np.where(lad_loc != fillvalues["tree_data"],(1.0 - lad_loc)*0.1,lad_loc)
     470   bad_loc = np.where(lad_loc != fillvalues["tree_data"],(1.0 - (lad_loc/(max(lad_loc.flatten())+0.01)))*0.1,lad_loc)
    462471   
    463472
     
    469478            if ( z[k] <= crown_center ):
    470479               r_test = np.sqrt( (x[i] - tree_location_x)**(2) + (y[j] - tree_location_y)**(2)  )
    471                if ( r_test <= radius ):
    472                   bad_loc[k,j,i] = 1.0
    473                if ( (r_test == 0.0) & (trunk_dia <= dx) ):
    474                    bad_loc[k,j,i] = radius**(2) * 3.14159265359
     480               if ( r_test == 0.0 ):
     481                  if ( trunk_dia <= dx ):
     482                     bad_loc[k,j,i] = radius**2 * 3.14159265359
     483                  else:
     484                     # WORKAROUND: divide remaining circle area over the 8 surrounding valid_pixels
     485                     bad_loc[k,j-1:j+2,i-1:i+2] = radius**2 * 3.14159265359 / 8.0
     486                     # for the central pixel fill the pixel
     487                     bad_loc[k,j,i] = dx**2
     488               #elif ( r_test <= radius ):
     489                  # TODO: calculate circle segment of grid points cut by the grid
     490 
    475491
    476492   return lad_loc, bad_loc, x, y, z, 0
    477493
    478494
    479 def process_patch(dz,patch_height,max_height_lad,patch_lai,alpha,beta):
     495def process_patch(dz,patch_height,vegetation_type,max_height_lad,patch_lai,alpha,beta):
    480496   
    481497#  Define fill values
    482498   fillvalues = {
    483499   "tree_data": float(-9999.0),
    484    "pch_index": int(-9999),
     500   "tree_id":   int(-9999),   
     501   "tree_type": np.byte(-127), 
     502   "pch_index": int(-9999)
    485503   }
    486504
     505   
    487506   phdz = patch_height[:,:] / dz
    488507   pch_index = np.where( (patch_height[:,:] != fillvalues["tree_data"]),phdz.astype(int)+1,int(-1))   
     
    520539                lad_loc[k,j,i] = 0.5 * ( pre_lad[k-1] + pre_lad[k] )
    521540
    522    return lad_loc, nz, 0
     541
     542   patch_id_2d = np.where(lad_loc[0,:,:] != fillvalues["tree_data"], 1, 0)
     543   patch_id = np.where(lad_loc != fillvalues["tree_data"], 1, 0)
     544   patch_type = np.where(lad_loc != fillvalues["tree_data"], np.byte(1), np.byte(-127))
     545
     546   num_ids = nd.label(patch_id_2d, output=patch_id_2d)
     547 
     548   for k in range(0,nz):
     549   
     550      vegetation_type_int = vegetation_type.astype(int)
     551      patch_id[k,:,:]   = np.where((patch_id_2d != 0) & (lad_loc[k,:,:] != fillvalues["tree_data"]), patch_id_2d, fillvalues["tree_id"])
     552      patch_type[k,:,:] = np.where((patch_id_2d != 0) & (lad_loc[k,:,:] != fillvalues["tree_data"]), np.byte(100+vegetation_type), fillvalues["tree_type"])
     553
     554   return lad_loc, patch_id, patch_type, nz, 0
    523555
    524556
     
    559591       self.bad_scale = bad_scale
    560592       self.dbh = dbh   
     593   
Note: See TracChangeset for help on using the changeset viewer.