Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

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

Last change on this file since 3648 was 3629, checked in by maronga, 6 years ago

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