source: palm/trunk/SCRIPTS/palm_csd_files/palm_csd_canopy_generator.py @ 3668

Last change on this file since 3668 was 3668, checked in by maronga, 3 years ago

removed most_methods circular and lookup. added improved version of palm_csd

  • Property svn:keywords set to Id
File size: 27.6 KB
Line 
1#!/usr/bin/env python3
2# -*- coding: utf-8 -*-
3#--------------------------------------------------------------------------------#
4# This file is part of the PALM model system.
5#
6# PALM is free software: you can redistribute it and/or modify it under the terms
7# of the GNU General Public License as published by the Free Software Foundation,
8# either version 3 of the License, or (at your option) any later version.
9#
10# PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12# A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13#
14# You should have received a copy of the GNU General Public License along with
15# PALM. If not, see <http://www.gnu.org/licenses/>.
16#
17# Copyright 1997-2018  Leibniz Universitaet Hannover
18#--------------------------------------------------------------------------------#
19#
20# Current revisions:
21# -----------------
22#
23#
24# Former revisions:
25# -----------------
26# $Id: palm_csd_canopy_generator.py 3668 2019-01-14 12:49:24Z maronga $
27# Various improvements and bugfixes
28#
29# 3629 2018-12-13 12:18:54Z maronga
30# Initial revision
31#
32# Description:
33# ------------
34# Canopy generator routines for creating 3D leaf and basal area densities for
35# single trees and tree patches
36#
37# @Author Bjoern Maronga (maronga@muk.uni-hannover.de)
38#------------------------------------------------------------------------------#
39import numpy as np
40import math
41import scipy.integrate as integrate
42
43def 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):
44
45
46#  Step 1: Create arrays for storing the data
47   max_canopy_height = max( max_tree_height,max_patch_height )
48 
49   zlad = np.arange(0,math.floor(max_canopy_height/dz)*dz+2*dz,dz)
50   zlad[1:] = zlad[1:] - 0.5 * dz
51
52   lad = np.ones((len(zlad),len(y),len(x)))
53   lad[:,:,:] = fill
54
55   bad = np.ones((len(zlad),len(y),len(x)))
56   bad[:,:,:] = fill
57   
58   ids = np.ones((len(zlad),len(y),len(x)))
59   ids[:,:,:] = fill
60   
61
62#  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   
63   number_of_trees_array = np.where( (tree_type.flatten() != fill) | (tree_dia.flatten() != fill) | (trunk_dia.flatten() != fill),1.0,fill)
64   number_of_trees = len( number_of_trees_array[number_of_trees_array == 1.0] )
65   dx = x[1] - x[0]
66       
67   valid_pixels = np.where( (tree_type[:,:] != fill) | (tree_dia[:,:] != fill) | (trunk_dia[:,:] != fill),True,False)
68   
69#  For each tree, create a small 3d array containing the LAD field for the individual tree
70   print("Start generating " +  str(number_of_trees) + " trees...")
71   tree_id_counter = 0
72   if (number_of_trees > 0):
73      for i in range(0,len(x)-1):
74         for j in range(0,len(y)-1):
75            if valid_pixels[j,i]:
76               tree_id_counter = tree_id_counter + 1
77#               print("   Processing tree No " +  str(tree_id_counter) + " ...", end="")
78               lad_loc, bad_loc, x_loc, y_loc, z_loc, status = process_single_tree(dx,dz,tree_type[j,i],tree_height[j,i],tree_lai[j,i],tree_dia[j,i],trunk_dia[j,i],season)
79               if ( np.any(lad_loc) != fill ):
80 
81#                 Calculate the position of the local 3d tree array within the full domain in order to achieve correct mapping and cutting off at the edges of the full domain
82                  lad_loc_nx = int(len(x_loc) / 2)
83                  lad_loc_ny = int(len(y_loc) / 2)
84                  lad_loc_nz = int(len(z_loc))
85                             
86                  odd_x = int(len(x_loc) % 2)
87                  odd_y = int(len(y_loc) % 2)
88               
89                  ind_l_x  = max(0,(i-lad_loc_nx))
90                  ind_l_y  = max(0,(j-lad_loc_ny))
91                  ind_r_x  = min(len(x)-1,i+lad_loc_nx-1+odd_x)
92                  ind_r_y  = min(len(y)-1,j+lad_loc_ny-1+odd_y)
93               
94                  out_l_x = ind_l_x - (i-lad_loc_nx)
95                  out_l_y = ind_l_y - (j-lad_loc_ny)
96                  out_r_x = len(x_loc)-1 + ind_r_x - (i+lad_loc_nx-1+odd_x)
97                  out_r_y = len(y_loc)-1 + ind_r_y - (j+lad_loc_ny-1+odd_y)               
98   
99
100                  lad[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,lad_loc[0:lad_loc_nz,out_l_y:out_r_y+1,out_l_x:out_r_x+1],lad[0:lad_loc_nz,ind_l_y:ind_r_y+1,ind_l_x:ind_r_x+1]) 
101                  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]) 
102                  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]) 
103
104
105#                  if ( status == 0 ):
106#                     status_char = " ok."
107#                  else:
108#                     status_char = " skipped."
109#                  print(status_char)
110
111               del lad_loc, x_loc, y_loc, z_loc, status
112   return lad, bad, ids, zlad
113
114
115def process_single_tree(dx,dz,tree_type,tree_height,tree_lai,tree_dia,trunk_dia,season):
116
117#  Set some parameters
118   sphere_extinction = 0.6
119   cone_extinction = 0.2
120   ml_n_low = 0.5
121   ml_n_high = 6.0
122   
123#  Populate look up table for tree species and their properties
124#  #1 Tree shapes were manually lookep up.
125#  #2 Crown h/w ratio - missing
126#  #3 Crown diameter based on Berlin tree statistics
127#  #4 Tree height based on Berlin tree statistics
128#  #5 Tree LAI summer - missing
129#  #6 Tree LAI winter - missing
130#  #7 Height of lad maximum - missing
131#  #8 Ratio LAD/BAD - missing
132#  #9 Trunk diameter at breast height from Berlin
133   default_trees = []
134                                                #1   #2   #3   #4    #5   #6   #7   #8     #9
135   default_trees.append(tree("Default",         1.0, 1.0, 4.0, 12.0, 3.0, 0.8, 0.6, 0.025, 0.35)) 
136   default_trees.append(tree("Abies",           3.0, 1.0, 4.0, 12.0, 3.0, 0.8, 0.6, 0.025, 0.80)) 
137   default_trees.append(tree("Acer",            1.0, 1.0, 7.0, 12.0, 3.0, 0.8, 0.6, 0.025, 0.80))
138   default_trees.append(tree("Aesculus",        1.0, 1.0, 7.0, 12.0, 3.0, 0.8, 0.6, 0.025, 1.00))
139   default_trees.append(tree("Ailanthus",       1.0, 1.0, 8.5, 13.5, 3.0, 0.8, 0.6, 0.025, 1.30))
140   default_trees.append(tree("Alnus",           3.0, 1.0, 6.0, 16.0, 3.0, 0.8, 0.6, 0.025, 1.20))
141   default_trees.append(tree("Amelanchier",     1.0, 1.0, 3.0,  4.0, 3.0, 0.8, 0.6, 0.025, 1.20))
142   default_trees.append(tree("Betula",          1.0, 1.0, 6.0, 14.0, 3.0, 0.8, 0.6, 0.025, 0.30))
143   default_trees.append(tree("Buxus",           1.0, 1.0, 4.0,  4.0, 3.0, 0.8, 0.6, 0.025, 0.90))
144   default_trees.append(tree("Calocedrus",      3.0, 1.0, 5.0, 10.0, 3.0, 0.8, 0.6, 0.025, 0.50))
145   default_trees.append(tree("Caragana",        1.0, 1.0, 3.5,  6.0, 3.0, 0.8, 0.6, 0.025, 0.90))
146   default_trees.append(tree("Carpinus",        1.0, 1.0, 6.0, 10.0, 3.0, 0.8, 0.6, 0.025, 0.70))
147   default_trees.append(tree("Carya",           1.0, 1.0, 5.0, 17.0, 3.0, 0.8, 0.6, 0.025, 0.80))
148   default_trees.append(tree("Castanea",        1.0, 1.0, 4.5,  7.0, 3.0, 0.8, 0.6, 0.025, 0.80))
149   default_trees.append(tree("Catalpa",         1.0, 1.0, 5.5,  6.5, 3.0, 0.8, 0.6, 0.025, 0.70))
150   default_trees.append(tree("Cedrus",          1.0, 1.0, 8.0, 13.0, 3.0, 0.8, 0.6, 0.025, 0.80))
151   default_trees.append(tree("Celtis",          1.0, 1.0, 6.0,  9.0, 3.0, 0.8, 0.6, 0.025, 0.80))
152   default_trees.append(tree("Cercidiphyllum",  1.0, 1.0, 3.0,  6.5, 3.0, 0.8, 0.6, 0.025, 0.80))
153   default_trees.append(tree("Cercis",          1.0, 1.0, 2.5,  7.5, 3.0, 0.8, 0.6, 0.025, 0.90))
154   default_trees.append(tree("Chamaecyparis",   5.0, 1.0, 3.5,  9.0, 3.0, 0.8, 0.6, 0.025, 0.70))
155   default_trees.append(tree("Cladrastis",      1.0, 1.0, 5.0, 10.0, 3.0, 0.8, 0.6, 0.025, 0.80))
156   default_trees.append(tree("Cornus",          1.0, 1.0, 4.5,  6.5, 3.0, 0.8, 0.6, 0.025, 1.20))
157   default_trees.append(tree("Corylus",         1.0, 1.0, 5.0,  9.0, 3.0, 0.8, 0.6, 0.025, 0.40))
158   default_trees.append(tree("Cotinus",         1.0, 1.0, 4.0,  4.0, 3.0, 0.8, 0.6, 0.025, 0.70))
159   default_trees.append(tree("Crataegus",       3.0, 1.0, 3.5,  6.0, 3.0, 0.8, 0.6, 0.025, 1.40))
160   default_trees.append(tree("Cryptomeria",     3.0, 1.0, 5.0, 10.0, 3.0, 0.8, 0.6, 0.025, 0.50))
161   default_trees.append(tree("Cupressocyparis", 3.0, 1.0, 3.0,  8.0, 3.0, 0.8, 0.6, 0.025, 0.40))
162   default_trees.append(tree("Cupressus",       3.0, 1.0, 5.0,  7.0, 3.0, 0.8, 0.6, 0.025, 0.40))
163   default_trees.append(tree("Cydonia",         1.0, 1.0, 2.0,  3.0, 3.0, 0.8, 0.6, 0.025, 0.90))
164   default_trees.append(tree("Davidia",         1.0, 1.0,10.0, 14.0, 3.0, 0.8, 0.6, 0.025, 0.40))
165   default_trees.append(tree("Elaeagnus",       1.0, 1.0, 6.5,  6.0, 3.0, 0.8, 0.6, 0.025, 1.20))
166   default_trees.append(tree("Euodia",          1.0, 1.0, 4.5,  6.0, 3.0, 0.8, 0.6, 0.025, 0.90))
167   default_trees.append(tree("Euonymus",        1.0, 1.0, 4.5,  6.0, 3.0, 0.8, 0.6, 0.025, 0.60))
168   default_trees.append(tree("Fagus",           1.0, 1.0,10.0, 12.5, 3.0, 0.8, 0.6, 0.025, 0.50))
169   default_trees.append(tree("Fraxinus",        1.0, 1.0, 5.5, 10.5, 3.0, 0.8, 0.6, 0.025, 1.60))
170   default_trees.append(tree("Ginkgo",          3.0, 1.0, 4.0,  8.5, 3.0, 0.8, 0.6, 0.025, 0.80))
171   default_trees.append(tree("Gleditsia",       1.0, 1.0, 6.5, 10.5, 3.0, 0.8, 0.6, 0.025, 0.60))
172   default_trees.append(tree("Gymnocladus",     1.0, 1.0, 5.5, 10.0, 3.0, 0.8, 0.6, 0.025, 0.80))
173   default_trees.append(tree("Hippophae",       1.0, 1.0, 9.5,  8.5, 3.0, 0.8, 0.6, 0.025, 0.80))
174   default_trees.append(tree("Ilex",            1.0, 1.0, 4.0,  7.5, 3.0, 0.8, 0.6, 0.025, 0.80))
175   default_trees.append(tree("Juglans",         1.0, 1.0, 7.0,  9.0, 3.0, 0.8, 0.6, 0.025, 0.50))
176   default_trees.append(tree("Juniperus",       5.0, 1.0, 3.0,  7.0, 3.0, 0.8, 0.6, 0.025, 0.90))
177   default_trees.append(tree("Koelreuteria",    1.0, 1.0, 3.5,  5.5, 3.0, 0.8, 0.6, 0.025, 0.50))
178   default_trees.append(tree("Laburnum",        1.0, 1.0, 3.0,  6.0, 3.0, 0.8, 0.6, 0.025, 0.60))
179   default_trees.append(tree("Larix",           3.0, 1.0, 7.0, 16.5, 3.0, 0.8, 0.6, 0.025, 0.60))
180   default_trees.append(tree("Ligustrum",       1.0, 1.0, 3.0,  6.0, 3.0, 0.8, 0.6, 0.025, 1.10))
181   default_trees.append(tree("Liquidambar",     3.0, 1.0, 3.0,  7.0, 3.0, 0.8, 0.6, 0.025, 0.30))
182   default_trees.append(tree("Liriodendron",    3.0, 1.0, 4.5,  9.5, 3.0, 0.8, 0.6, 0.025, 0.50))
183   default_trees.append(tree("Lonicera",        1.0, 1.0, 7.0,  9.0, 3.0, 0.8, 0.6, 0.025, 0.70))
184   default_trees.append(tree("Magnolia",        1.0, 1.0, 3.0,  5.0, 3.0, 0.8, 0.6, 0.025, 0.60))
185   default_trees.append(tree("Malus",           1.0, 1.0, 4.5,  5.0, 3.0, 0.8, 0.6, 0.025, 0.30))
186   default_trees.append(tree("Metasequoia",     5.0, 1.0, 4.5, 12.0, 3.0, 0.8, 0.6, 0.025, 0.50))
187   default_trees.append(tree("Morus",           1.0, 1.0, 7.5, 11.5, 3.0, 0.8, 0.6, 0.025, 1.00))
188   default_trees.append(tree("Ostrya",          1.0, 1.0, 2.0,  6.0, 3.0, 0.8, 0.6, 0.025, 1.00))
189   default_trees.append(tree("Parrotia",        1.0, 1.0, 7.0,  7.0, 3.0, 0.8, 0.6, 0.025, 0.30))
190   default_trees.append(tree("Paulownia",       1.0, 1.0, 4.0,  8.0, 3.0, 0.8, 0.6, 0.025, 0.40))
191   default_trees.append(tree("Phellodendron",   1.0, 1.0,13.5, 13.5, 3.0, 0.8, 0.6, 0.025, 0.50))
192   default_trees.append(tree("Picea",           3.0, 1.0, 3.0, 13.0, 3.0, 0.8, 0.6, 0.025, 0.90))
193   default_trees.append(tree("Pinus",           3.0, 1.0, 6.0, 16.0, 3.0, 0.8, 0.6, 0.025, 0.80))
194   default_trees.append(tree("Platanus",        1.0, 1.0,10.0, 14.5, 3.0, 0.8, 0.6, 0.025, 1.10))
195   default_trees.append(tree("Populus",         1.0, 1.0, 9.0, 20.0, 3.0, 0.8, 0.6, 0.025, 1.40))
196   default_trees.append(tree("Prunus",          1.0, 1.0, 5.0,  7.0, 3.0, 0.8, 0.6, 0.025, 1.60))
197   default_trees.append(tree("Pseudotsuga",     3.0, 1.0, 6.0, 17.5, 3.0, 0.8, 0.6, 0.025, 0.70))
198   default_trees.append(tree("Ptelea",          1.0, 1.0, 5.0,  4.0, 3.0, 0.8, 0.6, 0.025, 1.10))
199   default_trees.append(tree("Pterocaria",      1.0, 1.0,10.0, 12.0, 3.0, 0.8, 0.6, 0.025, 0.50))
200   default_trees.append(tree("Pterocarya",      1.0, 1.0,11.5, 14.5, 3.0, 0.8, 0.6, 0.025, 1.60))
201   default_trees.append(tree("Pyrus",           3.0, 1.0, 3.0,  6.0, 3.0, 0.8, 0.6, 0.025, 1.80))
202   default_trees.append(tree("Quercus",         1.0, 1.0, 8.0, 14.0, 3.1, 0.1, 0.6, 0.025, 0.40)) #
203   default_trees.append(tree("Rhamnus",         1.0, 1.0, 4.5,  4.5, 3.0, 0.8, 0.6, 0.025, 1.30))
204   default_trees.append(tree("Rhus",            1.0, 1.0, 7.0,  5.5, 3.0, 0.8, 0.6, 0.025, 0.50))
205   default_trees.append(tree("Robinia",         1.0, 1.0, 4.5, 13.5, 3.0, 0.8, 0.6, 0.025, 0.50))
206   default_trees.append(tree("Salix",           1.0, 1.0, 7.0, 14.0, 3.0, 0.8, 0.6, 0.025, 1.10))
207   default_trees.append(tree("Sambucus",        1.0, 1.0, 8.0,  6.0, 3.0, 0.8, 0.6, 0.025, 1.40))
208   default_trees.append(tree("Sasa",            1.0, 1.0,10.0, 25.0, 3.0, 0.8, 0.6, 0.025, 0.60))
209   default_trees.append(tree("Sequoiadendron",  5.0, 1.0, 5.5, 10.5, 3.0, 0.8, 0.6, 0.025, 1.60))
210   default_trees.append(tree("Sophora",         1.0, 1.0, 7.5, 10.0, 3.0, 0.8, 0.6, 0.025, 1.40))
211   default_trees.append(tree("Sorbus",          1.0, 1.0, 4.0,  7.0, 3.0, 0.8, 0.6, 0.025, 1.10))
212   default_trees.append(tree("Syringa",         1.0, 1.0, 4.5,  5.0, 3.0, 0.8, 0.6, 0.025, 0.60))
213   default_trees.append(tree("Tamarix",         1.0, 1.0, 6.0,  7.0, 3.0, 0.8, 0.6, 0.025, 0.50))
214   default_trees.append(tree("Taxodium",        5.0, 1.0, 6.0, 16.5, 3.0, 0.8, 0.6, 0.025, 0.60))
215   default_trees.append(tree("Taxus",           2.0, 1.0, 5.0,  7.5, 3.0, 0.8, 0.6, 0.025, 1.50))
216   default_trees.append(tree("Thuja",           3.0, 1.0, 3.5,  9.0, 3.0, 0.8, 0.6, 0.025, 0.70))
217   default_trees.append(tree("Tilia",           3.0, 1.0, 7.0, 12.5, 3.0, 0.8, 0.6, 0.025, 0.70))
218   default_trees.append(tree("Tsuga",           3.0, 1.0, 6.0, 10.5, 3.0, 0.8, 0.6, 0.025, 1.10))
219   default_trees.append(tree("Ulmus",           1.0, 1.0, 7.5, 14.0, 3.0, 0.8, 0.6, 0.025, 0.80))
220   default_trees.append(tree("Zelkova",         1.0, 1.0, 4.0,  5.5, 3.0, 0.8, 0.6, 0.025, 1.20))
221   default_trees.append(tree("Zenobia",         1.0, 1.0, 5.0,  5.0, 3.0, 0.8, 0.6, 0.025, 0.40))
222   
223#  Define fill values
224   fillvalues = {
225   "tree_data": float(-9999.0),
226   }
227   
228   
229#  Check for missing data in the input and set default values if needed
230   if ( tree_type == fillvalues["tree_data"] ):
231      tree_type = int(0)
232   else:
233      tree_type = int(tree_type)
234
235   if ( tree_height == fillvalues["tree_data"] ):
236      tree_height = default_trees[tree_type].height
237     
238   if ( tree_lai == fillvalues["tree_data"] ):
239      if (season == "summer"):
240         tree_lai = default_trees[tree_type].lai_summer
241
242      else:
243         tree_lai = default_trees[tree_type].lai_winter     
244
245   if ( tree_dia == fillvalues["tree_data"] ):
246      tree_dia = default_trees[tree_type].diameter
247     
248   if ( trunk_dia == fillvalues["tree_data"] ):
249      trunk_dia = default_trees[tree_type].dbh   
250 
251 
252#  Assign values that are not defined as user input from lookup table
253   tree_shape          = default_trees[tree_type].shape
254   tree_ratio          = default_trees[tree_type].ratio
255   lad_max_height      = default_trees[tree_type].lad_max_height
256   bad_scale           = default_trees[tree_type].bad_scale
257
258   #print("Tree input parameters:")
259   #print("----------------------")
260   #print("type:           " + str(default_trees[tree_type].species) )
261   #print("height:         " + str(tree_height))
262   #print("lai:            " + str(tree_lai))
263   #print("crown diameter: " + str(tree_dia))
264   #print("trunk diameter: " + str(trunk_dia))
265   #print("shape: " + str(tree_shape))
266   #print("height/width: " + str(tree_ratio))
267
268#  Calculate crown height and height of the crown center
269   crown_height   = tree_ratio * tree_dia
270   if ( crown_height > tree_height ):
271      crown_height = tree_height
272
273   crown_center   = tree_height - crown_height * 0.5
274
275
276#  Calculate height of maximum LAD
277   z_lad_max = lad_max_height * tree_height
278   
279#  Calculate the maximum LAD after Lalic and Mihailovic (2004)
280   lad_max_part_1 = integrate.quad(lambda z: ( ( tree_height - z_lad_max ) / ( tree_height - z ) ) ** (ml_n_high) * np.exp( ml_n_high * (1.0 - ( tree_height - z_lad_max ) / ( tree_height - z ) ) ), 0.0, z_lad_max)
281   lad_max_part_2 = integrate.quad(lambda z: ( ( tree_height - z_lad_max ) / ( tree_height - z ) ) ** (ml_n_low) * np.exp( ml_n_low * (1.0 - ( tree_height - z_lad_max ) / ( tree_height - z ) ) ), z_lad_max, tree_height)
282   
283   lad_max = tree_lai / (lad_max_part_1[0] + lad_max_part_2[0])
284   
285   
286#  Define position of tree and its output domain
287   nx = int(tree_dia / dx) + 2
288   nz = int(tree_height / dz) + 2
289   
290#  Add one grid point if diameter is an odd value   
291   if ( (tree_dia % 2.0) != 0.0 ):
292      nx = nx + 1
293
294
295#  Create local domain of the tree's LAD
296   x = np.arange(0,nx*dx,dx)
297   x[:] = x[:] - 0.5 * dx
298   y = x
299
300   z = np.arange(0,nz*dz,dz)
301   z[1:] = z[1:] - 0.5 * dz
302
303#  Define center of the tree position inside the local LAD domain
304   tree_location_x = x[nx/2]
305   tree_location_y = y[nx/2]
306
307
308#  Calculate LAD profile after Lalic and Mihailovic (2004). Will be later used for normalization
309   lad_profile     =  np.arange(0,nz,1.0)
310   lad_profile[:]  = 0.0
311
312   for k in range(1,nz-1):
313       if ( (z[k] > 0.0) & (z[k] < z_lad_max) ):
314          n = ml_n_high
315       else:
316          n = ml_n_low
317     
318       lad_profile[k] =  lad_max * ( ( tree_height - z_lad_max ) / ( tree_height - z[k] ) ) ** (n) * np.exp( n * (1.0 - ( tree_height - z_lad_max ) / ( tree_height - z[k] ) ) )
319
320#  Create lad array and populate according to the specific tree shape. This is still experimental
321   lad_loc = np.ones((nz,nx,nx))
322   lad_loc[:,:,:] = fillvalues["tree_data"]
323   bad_loc = np.copy(lad_loc)
324
325
326#  For very small trees, no LAD is calculated
327   if ( tree_height <=  (0.5*dz) ):
328        print("    Shallow tree found. Action: ignore.")
329        return lad_loc, bad_loc, x, y, z, 1
330   
331   
332#  Branch for spheres and ellipsoids. A symmetric LAD sphere is created assuming an LAD extinction towards the center of the tree, representing the effect of sunlight extinction and thus
333#  less leaf mass inside the tree crown. Extinction coefficients are experimental.
334   if ( tree_shape == 1 ):
335      for i in range(0,nx-1):
336         for j in range(0,nx-1):
337            for k in range(0,nz):
338               r_test = np.sqrt( (x[i] - tree_location_x)**(2)/(tree_dia * 0.5)**(2) + (y[j] - tree_location_y)**(2)/(tree_dia * 0.5)**2 + (z[k] - crown_center)**(2)/(crown_height * 0.5)**(2))
339               if ( r_test <= 1.0 ):
340                  lad_loc[k,j,i] = lad_max * np.exp( - sphere_extinction * (1.0 - r_test) ) 
341               else:
342                  lad_loc[k,j,i] = fillvalues["tree_data"]
343
344            if ( np.any( lad_loc[:,j,i] != fillvalues["tree_data"]) ):
345               lad_loc[0,j,i] = 0.0
346
347
348#  Branch for cylinder shapes
349   if ( tree_shape == 2 ):
350      k_min = int((crown_center - crown_height * 0.5) / dz)
351      k_max = int((crown_center + crown_height * 0.5) / dz)
352      for i in range(0,nx-1):
353         for j in range(0,nx-1):
354            for k in range(k_min,k_max):
355               r_test = np.sqrt( (x[i] - tree_location_x)**(2)/(tree_dia * 0.5)**(2) + (y[j] - tree_location_y)**(2)/(tree_dia * 0.5)**(2))
356               if ( r_test <= 1.0 ):
357                  r_test3 = np.sqrt( (z[k] - crown_center)**(2)/(crown_height * 0.5)**(2))
358                  lad_loc[k,j,i] = lad_max * np.exp ( - sphere_extinction * (1.0 - max(r_test,r_test3) ) )                 
359               else:
360                  lad_loc[k,j,i] = fillvalues["tree_data"]
361
362            if ( np.any( lad_loc[:,j,i] != fillvalues["tree_data"]) ):
363               lad_loc[0,j,i] = 0.0
364
365#  Branch for cone shapes
366   if ( tree_shape == 3 ):
367      k_min = int((crown_center - crown_height * 0.5) / dz)
368      k_max = int((crown_center + crown_height * 0.5) / dz)
369      for i in range(0,nx-1):
370         for j in range(0,nx-1):
371            for k in range(k_min,k_max):
372               k_rel = k - k_min
373               r_test = (x[i] - tree_location_x)**(2) + (y[j] - tree_location_y)**(2) - ( (tree_dia * 0.5)**(2) / crown_height**(2) ) * ( z[k_rel] - crown_height)**(2)
374               if ( r_test <= 0.0 ):
375                  r_test2 = np.sqrt( (x[i] - tree_location_x)**(2)/(tree_dia * 0.5)**(2) + (y[j] - tree_location_y)**(2)/(tree_dia * 0.5)**(2))
376                  r_test3 = np.sqrt( (z[k] - crown_center)**(2)/(crown_height * 0.5)**(2))
377                  lad_loc[k,j,i] = lad_max * np.exp ( - cone_extinction * (1.0 - max((r_test+1.0),r_test2,r_test3)) )         
378               else:
379                  lad_loc[k,j,i] = fillvalues["tree_data"]
380
381            if ( np.any( lad_loc[:,j,i] != fillvalues["tree_data"]) ):
382               lad_loc[0,j,i] = 0.0
383 
384#  Branch for inverted cone shapes. TODO: what is r_test2 and r_test3 used for? Debugging needed!
385   if ( tree_shape == 4 ):
386      k_min = int((crown_center - crown_height * 0.5) / dz)
387      k_max = int((crown_center + crown_height * 0.5) / dz)
388      for i in range(0,nx-1):
389         for j in range(0,nx-1):
390            for k in range(k_min,k_max):
391               k_rel = k_max - k
392               r_test = (x[i] - tree_location_x)**(2) + (y[j] - tree_location_y)**(2) - ( (tree_dia * 0.5)**(2) / crown_height**(2) ) * ( z[k_rel] - crown_height)**(2)
393               if ( r_test <= 0.0 ):
394                  r_test2 = np.sqrt( (x[i] - tree_location_x)**(2)/(tree_dia * 0.5)**(2) + (y[j] - tree_location_y)**(2)/(tree_dia * 0.5)**(2))
395                  r_test3 = np.sqrt( (z[k] - crown_center)**(2)/(crown_height * 0.5)**(2))
396                  lad_loc[k,j,i] = lad_max * np.exp ( - cone_extinction * ( - r_test ) )         
397               else:
398                  lad_loc[k,j,i] = fillvalues["tree_data"]
399
400            if ( np.any( lad_loc[:,j,i] != fillvalues["tree_data"]) ):
401               lad_loc[0,j,i] = 0.0
402 
403 
404#  Branch for paraboloid shapes
405   if ( tree_shape == 5 ):
406      k_min = int((crown_center - crown_height * 0.5) / dz)
407      k_max = int((crown_center + crown_height * 0.5) / dz)
408      for i in range(0,nx-1):
409         for j in range(0,nx-1):
410            for k in range(k_min,k_max):
411               k_rel = k - k_min
412               r_test = ((x[i] - tree_location_x)**(2) + (y[j] - tree_location_y)**(2)) * crown_height / (tree_dia * 0.5)**(2) - z[k_rel]
413               if ( r_test <= 0.0 ):
414                   lad_loc[k,j,i] = lad_max * np.exp ( - cone_extinction * (- r_test) )               
415               else:
416                  lad_loc[k,j,i] = fillvalues["tree_data"]
417
418            if ( np.any( lad_loc[:,j,i] != fillvalues["tree_data"]) ):
419               lad_loc[0,j,i] = 0.0
420 
421 
422 
423#  Branch for inverted paraboloid shapes
424   if ( tree_shape == 6 ):
425      k_min = int((crown_center - crown_height * 0.5) / dz)
426      k_max = int((crown_center + crown_height * 0.5) / dz)
427      for i in range(0,nx-1):
428         for j in range(0,nx-1):
429            for k in range(k_min,k_max):
430               k_rel = k_max - k
431               r_test = ((x[i] - tree_location_x)**(2) + (y[j] - tree_location_y)**(2)) * crown_height / (tree_dia * 0.5)**(2) - z[k_rel]
432               if ( r_test <= 0.0 ):
433                   lad_loc[k,j,i] = lad_max * np.exp ( - cone_extinction * (- r_test) )               
434               else:
435                  lad_loc[k,j,i] = fillvalues["tree_data"]
436
437            if ( np.any( lad_loc[:,j,i] != fillvalues["tree_data"]) ):
438               lad_loc[0,j,i] = 0.0 
439 
440 
441#  Normalize the LAD profile so that the vertically integrated Lalic and Mihailovic (2004) is reproduced by the LAD array. Deactivated for now.
442   #for i in range(0,nx-1):
443      #for j in range(0,nx-1):
444         #lad_clean = np.where(lad_loc[:,j,i] == fillvalues["tree_data"],0.0,lad_loc[:,j,i])
445         #lai_from_int = integrate.simps (lad_clean, z)
446         #print(lai_from_int)
447         #for k in range(0,nz):
448            #if ( np.any(lad_loc[k,j,i] > 0.0) ): 
449               #lad_loc[k,j,i] = np.where((lad_loc[k,j,i] != fillvalues["tree_data"]),lad_loc[k,j,i] / lai_from_int * lad_profile[k],lad_loc[k,j,i])
450
451
452#  Dreate BAD array and populate
453   bad_loc = np.where(lad_loc != fillvalues["tree_data"],lad_loc*0.01,lad_loc)
454   
455
456#  Overwrite grid cells that are occupied by the tree trunk
457   radius = trunk_dia * 0.5
458   for i in range(0,nx-1):
459      for j in range(0,nx-1):
460         for k in range(0,nz):
461            if ( z[k] <= crown_center ):
462               r_test = np.sqrt( (x[i] - tree_location_x)**(2) + (y[j] - tree_location_y)**(2)  )
463               if ( r_test <= radius ):
464                  bad_loc[k,j,i] = 1.0
465               if ( (r_test == 0.0) & (trunk_dia <= dx) ):
466                   bad_loc[k,j,i] = radius**(2) * 3.14159265359
467
468   return lad_loc, bad_loc, x, y, z, 0
469
470
471def process_patch(dz,patch_height,max_height_lad,patch_lai,alpha,beta):
472   
473#  Define fill values
474   fillvalues = {
475   "tree_data": float(-9999.0),
476   "pch_index": int(-9999),
477   }
478
479   phdz = patch_height[:,:] / dz
480   pch_index = np.where( (patch_height[:,:] != fillvalues["tree_data"]),phdz.astype(int)+1,int(-1))   
481   pch_index = np.where( (pch_index[:,:] == 0) ,fillvalues["pch_index"],pch_index[:,:])
482   pch_index = np.where( (pch_index[:,:] == -1) ,0,pch_index[:,:])
483
484   max_canopy_height = max(max(patch_height.flatten()),max_height_lad)
485
486   z = np.arange(0,math.floor(max_canopy_height/dz)*dz+2*dz,dz)
487   
488   z[1:] = z[1:] - 0.5 * dz
489
490   nz = len(z)
491   ny = len(patch_height[:,0])
492   nx = len(patch_height[0,:])
493
494   pre_lad = np.ones((nz))
495   pre_lad[:] = 0.0
496   lad_loc = np.ones( (nz,ny,nx) )
497   lad_loc[:,:,:] = fillvalues["tree_data"]
498   
499   for i in range(0,nx-1):
500      for j in range(0,ny-1):
501          int_bpdf = 0.0
502          if ( (patch_height[j,i] != fillvalues["tree_data"]) & (patch_height[j,i] >= (0.5*dz)) ):
503             for k in range(1,pch_index[j,i]):
504                int_bpdf = int_bpdf + ( ( ( z[k] / patch_height[j,i] )**( alpha - 1 ) ) * ( ( 1.0 - ( z[k] / patch_height[j,i] ) )**(beta - 1 ) ) * ( dz / patch_height[j,i] ) )
505
506             for k in range(1,pch_index[j,i]):
507                pre_lad[k] =  patch_lai[j,i] * ( ( ( dz*k / patch_height[j,i] )**( alpha - 1.0 ) ) * ( ( 1.0 - ( dz*k / patch_height[j,i] ) )**(beta - 1.0 ) ) / int_bpdf ) / patch_height[j,i]
508
509             lad_loc[0,j,i] = pre_lad[0]
510             
511             for k in range(0,pch_index[j,i]):           
512                lad_loc[k,j,i] = 0.5 * ( pre_lad[k-1] + pre_lad[k] )
513
514   return lad_loc, nz, 0
515
516
517# CLASS TREE
518#
519# Default tree geometrical parameters:
520#
521# species: name of the tree type
522#
523# shape: defines the general shape of the tree and can be one of the following types:
524# 1.0 sphere or ellipsoid
525# 2.0 cylinder
526# 3.0 cone
527# 4.0 inverted cone
528# 5.0 paraboloid (rounded cone)
529# 6.0 inverted paraboloid (invertes rounded cone)
530
531# ratio:  ratio of maximum crown height to the maximum crown diameter
532# diameter: default crown diameter (m)
533# height:   default total height of the tree including trunk (m)
534# lai_summer: default leaf area index fully leafed
535# lai_winter: default winter-teim leaf area index
536# lad_max: default maximum leaf area density (m2/m3)
537# lad_max_height: default height where the leaf area density is maximum relative to total tree height
538# bad_scale: ratio of basal area in the crown area to the leaf area
539# dbh: default trunk diameter at breast height (1.4 m) (m)
540#
541class tree:
542    def __init__(self, species, shape, ratio, diameter, height, lai_summer, lai_winter, lad_max_height, bad_scale, dbh): 
543       self.species = species
544       self.shape = shape
545       self.ratio = ratio
546       self.diameter = diameter
547       self.height = height
548       self.lai_summer = lai_summer
549       self.lai_winter = lai_winter
550       self.lad_max_height = lad_max_height
551       self.bad_scale = bad_scale
552       self.dbh = dbh   
Note: See TracBrowser for help on using the repository browser.