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

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

fix for last commit
File size: 27.7 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# 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="")
81               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)
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
118def process_single_tree(dx,dz,tree_type,tree_height,tree_lai,tree_dia,trunk_dia,season):
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
238   if ( tree_height == fillvalues["tree_data"] ):
239      tree_height = default_trees[tree_type].height
240     
241   if ( tree_lai == fillvalues["tree_data"] ):
242      if (season == "summer"):
243         tree_lai = default_trees[tree_type].lai_summer
244
245      else:
246         tree_lai = default_trees[tree_type].lai_winter     
247
248   if ( tree_dia == fillvalues["tree_data"] ):
249      tree_dia = default_trees[tree_type].diameter
250     
251   if ( trunk_dia == fillvalues["tree_data"] ):
252      trunk_dia = default_trees[tree_type].dbh   
253 
254 
255#  Assign values that are not defined as user input from lookup table
256   tree_shape          = default_trees[tree_type].shape
257   tree_ratio          = default_trees[tree_type].ratio
258   lad_max_height      = default_trees[tree_type].lad_max_height
259   bad_scale           = default_trees[tree_type].bad_scale
260
261   #print("Tree input parameters:")
262   #print("----------------------")
263   #print("type:           " + str(default_trees[tree_type].species) )
264   #print("height:         " + str(tree_height))
265   #print("lai:            " + str(tree_lai))
266   #print("crown diameter: " + str(tree_dia))
267   #print("trunk diameter: " + str(trunk_dia))
268   #print("shape: " + str(tree_shape))
269   #print("height/width: " + str(tree_ratio))
270
271#  Calculate crown height and height of the crown center
272   crown_height   = tree_ratio * tree_dia
273   if ( crown_height > tree_height ):
274      crown_height = tree_height
275
276   crown_center   = tree_height - crown_height * 0.5
277
278
279#  Calculate height of maximum LAD
280   z_lad_max = lad_max_height * tree_height
281   
282#  Calculate the maximum LAD after Lalic and Mihailovic (2004)
283   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)
284   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)
285   
286   lad_max = tree_lai / (lad_max_part_1[0] + lad_max_part_2[0])
287   
288   
289#  Define position of tree and its output domain
290   nx = int(tree_dia / dx) + 2
291   nz = int(tree_height / dz) + 2
292   
293#  Add one grid point if diameter is an odd value   
294   if ( (tree_dia % 2.0) != 0.0 ):
295      nx = nx + 1
296
297
298#  Create local domain of the tree's LAD
299   x = np.arange(0,nx*dx,dx)
300   x[:] = x[:] - 0.5 * dx
301   y = x
302
303   z = np.arange(0,nz*dz,dz)
304   z[1:] = z[1:] - 0.5 * dz
305
306#  Define center of the tree position inside the local LAD domain
307   tree_location_x = x[int(nx/2)]
308   tree_location_y = y[int(nx/2)]
309
310
311#  Calculate LAD profile after Lalic and Mihailovic (2004). Will be later used for normalization
312   lad_profile     =  np.arange(0,nz,1.0)
313   lad_profile[:]  = 0.0
314
315   for k in range(1,nz-1):
316       if ( (z[k] > 0.0) & (z[k] < z_lad_max) ):
317          n = ml_n_high
318       else:
319          n = ml_n_low
320     
321       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] ) ) )
322
323#  Create lad array and populate according to the specific tree shape. This is still experimental
324   lad_loc = np.ones((nz,nx,nx))
325   lad_loc[:,:,:] = fillvalues["tree_data"]
326   bad_loc = np.copy(lad_loc)
327
328
329#  For very small trees, no LAD is calculated
330   if ( tree_height <=  (0.5*dz) ):
331        print("    Shallow tree found. Action: ignore.")
332        return lad_loc, bad_loc, x, y, z, 1
333   
334   
335#  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
336#  less leaf mass inside the tree crown. Extinction coefficients are experimental.
337   if ( tree_shape == 1 ):
338      for i in range(0,nx-1):
339         for j in range(0,nx-1):
340            for k in range(0,nz):
341               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))
342               if ( r_test <= 1.0 ):
343                  lad_loc[k,j,i] = lad_max * np.exp( - sphere_extinction * (1.0 - r_test) ) 
344               else:
345                  lad_loc[k,j,i] = fillvalues["tree_data"]
346
347            if ( np.any( lad_loc[:,j,i] != fillvalues["tree_data"]) ):
348               lad_loc[0,j,i] = 0.0
349
350
351#  Branch for cylinder shapes
352   if ( tree_shape == 2 ):
353      k_min = int((crown_center - crown_height * 0.5) / dz)
354      k_max = int((crown_center + crown_height * 0.5) / dz)
355      for i in range(0,nx-1):
356         for j in range(0,nx-1):
357            for k in range(k_min,k_max):
358               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))
359               if ( r_test <= 1.0 ):
360                  r_test3 = np.sqrt( (z[k] - crown_center)**(2)/(crown_height * 0.5)**(2))
361                  lad_loc[k,j,i] = lad_max * np.exp ( - sphere_extinction * (1.0 - max(r_test,r_test3) ) )                 
362               else:
363                  lad_loc[k,j,i] = fillvalues["tree_data"]
364
365            if ( np.any( lad_loc[:,j,i] != fillvalues["tree_data"]) ):
366               lad_loc[0,j,i] = 0.0
367
368#  Branch for cone shapes
369   if ( tree_shape == 3 ):
370      k_min = int((crown_center - crown_height * 0.5) / dz)
371      k_max = int((crown_center + crown_height * 0.5) / dz)
372      for i in range(0,nx-1):
373         for j in range(0,nx-1):
374            for k in range(k_min,k_max):
375               k_rel = k - k_min
376               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)
377               if ( r_test <= 0.0 ):
378                  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))
379                  r_test3 = np.sqrt( (z[k] - crown_center)**(2)/(crown_height * 0.5)**(2))
380                  lad_loc[k,j,i] = lad_max * np.exp ( - cone_extinction * (1.0 - max((r_test+1.0),r_test2,r_test3)) )         
381               else:
382                  lad_loc[k,j,i] = fillvalues["tree_data"]
383
384            if ( np.any( lad_loc[:,j,i] != fillvalues["tree_data"]) ):
385               lad_loc[0,j,i] = 0.0
386 
387#  Branch for inverted cone shapes. TODO: what is r_test2 and r_test3 used for? Debugging needed!
388   if ( tree_shape == 4 ):
389      k_min = int((crown_center - crown_height * 0.5) / dz)
390      k_max = int((crown_center + crown_height * 0.5) / dz)
391      for i in range(0,nx-1):
392         for j in range(0,nx-1):
393            for k in range(k_min,k_max):
394               k_rel = k_max - k
395               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)
396               if ( r_test <= 0.0 ):
397                  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))
398                  r_test3 = np.sqrt( (z[k] - crown_center)**(2)/(crown_height * 0.5)**(2))
399                  lad_loc[k,j,i] = lad_max * np.exp ( - cone_extinction * ( - r_test ) )         
400               else:
401                  lad_loc[k,j,i] = fillvalues["tree_data"]
402
403            if ( np.any( lad_loc[:,j,i] != fillvalues["tree_data"]) ):
404               lad_loc[0,j,i] = 0.0
405 
406 
407#  Branch for paraboloid shapes
408   if ( tree_shape == 5 ):
409      k_min = int((crown_center - crown_height * 0.5) / dz)
410      k_max = int((crown_center + crown_height * 0.5) / dz)
411      for i in range(0,nx-1):
412         for j in range(0,nx-1):
413            for k in range(k_min,k_max):
414               k_rel = k - k_min
415               r_test = ((x[i] - tree_location_x)**(2) + (y[j] - tree_location_y)**(2)) * crown_height / (tree_dia * 0.5)**(2) - z[k_rel]
416               if ( r_test <= 0.0 ):
417                   lad_loc[k,j,i] = lad_max * np.exp ( - cone_extinction * (- r_test) )               
418               else:
419                  lad_loc[k,j,i] = fillvalues["tree_data"]
420
421            if ( np.any( lad_loc[:,j,i] != fillvalues["tree_data"]) ):
422               lad_loc[0,j,i] = 0.0
423 
424 
425 
426#  Branch for inverted paraboloid shapes
427   if ( tree_shape == 6 ):
428      k_min = int((crown_center - crown_height * 0.5) / dz)
429      k_max = int((crown_center + crown_height * 0.5) / dz)
430      for i in range(0,nx-1):
431         for j in range(0,nx-1):
432            for k in range(k_min,k_max):
433               k_rel = k_max - k
434               r_test = ((x[i] - tree_location_x)**(2) + (y[j] - tree_location_y)**(2)) * crown_height / (tree_dia * 0.5)**(2) - z[k_rel]
435               if ( r_test <= 0.0 ):
436                   lad_loc[k,j,i] = lad_max * np.exp ( - cone_extinction * (- r_test) )               
437               else:
438                  lad_loc[k,j,i] = fillvalues["tree_data"]
439
440            if ( np.any( lad_loc[:,j,i] != fillvalues["tree_data"]) ):
441               lad_loc[0,j,i] = 0.0 
442 
443 
444#  Normalize the LAD profile so that the vertically integrated Lalic and Mihailovic (2004) is reproduced by the LAD array. Deactivated for now.
445   #for i in range(0,nx-1):
446      #for j in range(0,nx-1):
447         #lad_clean = np.where(lad_loc[:,j,i] == fillvalues["tree_data"],0.0,lad_loc[:,j,i])
448         #lai_from_int = integrate.simps (lad_clean, z)
449         #print(lai_from_int)
450         #for k in range(0,nz):
451            #if ( np.any(lad_loc[k,j,i] > 0.0) ): 
452               #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])
453
454
455#  Dreate BAD array and populate
456   bad_loc = np.where(lad_loc != fillvalues["tree_data"],lad_loc*0.01,lad_loc)
457   
458
459#  Overwrite grid cells that are occupied by the tree trunk
460   radius = trunk_dia * 0.5
461   for i in range(0,nx-1):
462      for j in range(0,nx-1):
463         for k in range(0,nz):
464            if ( z[k] <= crown_center ):
465               r_test = np.sqrt( (x[i] - tree_location_x)**(2) + (y[j] - tree_location_y)**(2)  )
466               if ( r_test <= radius ):
467                  bad_loc[k,j,i] = 1.0
468               if ( (r_test == 0.0) & (trunk_dia <= dx) ):
469                   bad_loc[k,j,i] = radius**(2) * 3.14159265359
470
471   return lad_loc, bad_loc, x, y, z, 0
472
473
474def process_patch(dz,patch_height,max_height_lad,patch_lai,alpha,beta):
475   
476#  Define fill values
477   fillvalues = {
478   "tree_data": float(-9999.0),
479   "pch_index": int(-9999),
480   }
481
482   phdz = patch_height[:,:] / dz
483   pch_index = np.where( (patch_height[:,:] != fillvalues["tree_data"]),phdz.astype(int)+1,int(-1))   
484   pch_index = np.where( (pch_index[:,:] == 0) ,fillvalues["pch_index"],pch_index[:,:])
485   pch_index = np.where( (pch_index[:,:] == -1) ,0,pch_index[:,:])
486
487   max_canopy_height = max(max(patch_height.flatten()),max_height_lad)
488
489   z = np.arange(0,math.floor(max_canopy_height/dz)*dz+2*dz,dz)
490   
491   z[1:] = z[1:] - 0.5 * dz
492
493   nz = len(z)
494   ny = len(patch_height[:,0])
495   nx = len(patch_height[0,:])
496
497   pre_lad = np.ones((nz))
498   pre_lad[:] = 0.0
499   lad_loc = np.ones( (nz,ny,nx) )
500   lad_loc[:,:,:] = fillvalues["tree_data"]
501   
502   for i in range(0,nx-1):
503      for j in range(0,ny-1):
504          int_bpdf = 0.0
505          if ( (patch_height[j,i] != fillvalues["tree_data"]) & (patch_height[j,i] >= (0.5*dz)) ):
506             for k in range(1,pch_index[j,i]):
507                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] ) )
508
509             for k in range(1,pch_index[j,i]):
510                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]
511
512             lad_loc[0,j,i] = pre_lad[0]
513             
514             for k in range(0,pch_index[j,i]):           
515                lad_loc[k,j,i] = 0.5 * ( pre_lad[k-1] + pre_lad[k] )
516
517   return lad_loc, nz, 0
518
519
520# CLASS TREE
521#
522# Default tree geometrical parameters:
523#
524# species: name of the tree type
525#
526# shape: defines the general shape of the tree and can be one of the following types:
527# 1.0 sphere or ellipsoid
528# 2.0 cylinder
529# 3.0 cone
530# 4.0 inverted cone
531# 5.0 paraboloid (rounded cone)
532# 6.0 inverted paraboloid (invertes rounded cone)
533# 
534# ratio:  ratio of maximum crown height to the maximum crown diameter
535# diameter: default crown diameter (m)
536# height:   default total height of the tree including trunk (m)
537# lai_summer: default leaf area index fully leafed
538# lai_winter: default winter-teim leaf area index
539# lad_max: default maximum leaf area density (m2/m3)
540# lad_max_height: default height where the leaf area density is maximum relative to total tree height
541# bad_scale: ratio of basal area in the crown area to the leaf area
542# dbh: default trunk diameter at breast height (1.4 m) (m)
543#
544class tree:
545    def __init__(self, species, shape, ratio, diameter, height, lai_summer, lai_winter, lad_max_height, bad_scale, dbh): 
546       self.species = species
547       self.shape = shape
548       self.ratio = ratio
549       self.diameter = diameter
550       self.height = height
551       self.lai_summer = lai_summer
552       self.lai_winter = lai_winter
553       self.lad_max_height = lad_max_height
554       self.bad_scale = bad_scale
555       self.dbh = dbh   
Note: See TracBrowser for help on using the repository browser.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |