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 @ 4250

Last change on this file since 4250 was 4023, checked in by maronga, 5 years ago

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