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

Last change on this file was 4843, checked in by raasch, 21 months ago

local namelist parameter added to switch off the module although the respective module namelist appears in the namelist file, further copyright updates

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