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

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

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