Home


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

source: palm/trunk/SCRIPTS/palm_csd @ 3677

Last change on this file since 3677 was 3668, checked in by maronga, 6 years ago

removed most_methods circular and lookup. added improved version of palm_csd
  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 54.6 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 3668 2019-01-14 12:49:24Z moh.hefny $
27# Various improvements and bugfixes
28#
29# 3629 2018-12-13 12:18:54Z maronga
30# Added canopy generator calls. Some improvements
31#
32# 3567 2018-11-27 13:59:21Z maronga
33# Initial revisions
34#
35# Description:
36# ------------
37# Processing tool for creating PIDS conform static drivers from rastered NetCDF
38# input
39#
40# @Author Bjoern Maronga (maronga@muk.uni-hannover.de)
41#
42# @todo Make input files optional
43# @todo Allow for ASCII input of terrain height and building height
44# @todo Modularize reading config file
45#------------------------------------------------------------------------------#
46
47from palm_csd_files.palm_csd_netcdf_interface import *
48from palm_csd_files.palm_csd_tools import *
49from palm_csd_files.palm_csd_canopy_generator import *
50import numpy as np
51
52   
53def read_config_file():
54   
55   import configparser
56   from math import floor
57   import numpy as np
58   import os
59   import subprocess as sub
60   import sys
61
62#  Check if configuration files exists and quit otherwise
63   input_config = ".csd.config"
64   for i in range(1,len(sys.argv)): 
65      input_config = str(sys.argv[i])
66       
67   config = configparser.RawConfigParser(allow_no_value=True)
68
69   if ( os.path.isfile(input_config) == False ):
70      print ("Error. No configuration file " + input_config + " found.")
71      raise SystemExit   
72
73   config.read(input_config)
74
75
76#  Definition of settings
77
78   global settings_bridge_width
79   global settings_lai_roof_intensive
80   global settings_lai_roof_extensive
81   global settings_season
82   global settings_lai_low_default
83   global settings_lai_high_default
84   global settings_patch_height_default
85   global settings_lai_alpha
86   global settings_lai_beta
87   global ndomains
88
89#  Definition of global configuration parameters
90   global global_acronym
91   global global_angle
92   global global_author
93   global global_campaign
94   global global_comment
95   global global_contact
96   global global_data_content 
97   global global_dependencies
98   global global_institution
99   global global_keywords
100   global global_location
101   global global_palm_version
102   global global_references
103   global global_site
104   global global_source
105
106   global  path_out
107   global filename_out
108   global version_out
109
110
111#  Definition of domain parameters
112   global domain_names
113   global domain_px
114   global domain_x0
115   global domain_y0
116   global domain_x1
117   global domain_y1
118   global domain_nx
119   global domain_ny
120   global domain_dz
121   global domain_3d
122   global domain_high_vegetation
123   global domain_ip
124   global domain_za
125   global domain_parent
126   global domain_green_roofs
127   global domain_street_trees
128   global domain_canopy_patches
129
130#  Definition of input data parameters
131   global input_names
132   global input_px
133 
134       
135   global input_file_x
136   global input_file_y
137   global input_file_x_UTM
138   global input_file_y_UTM   
139   global input_file_lat
140   global input_file_lon 
141   global input_file_zt
142   global input_file_buildings_2d
143   global input_file_bridges_2d
144   global input_file_building_id
145   global input_file_bridges_id   
146   global input_file_building_type
147   global input_file_building_type
148   global input_file_lai 
149   global input_file_vegetation_type
150   global input_file_vegetation_height
151   global input_file_pavement_type
152   global input_file_water_type
153   global input_file_street_type
154   global input_file_street_crossings
155   global input_file_soil_type
156   global input_file_vegetation_on_roofs
157   global input_file_tree_crown_diameter
158   global input_file_tree_height
159   global input_file_tree_trunk_diameter
160   global input_file_tree_type
161   global input_file_patch_height
162
163   global zt_all
164   global zt_min
165
166   settings_bridge_width = 3.0
167   settings_lai_roof_intensive = 0.5
168   settings_lai_roof_extensive = 1.0
169   settings_season = "summer"
170   settings_lai_high_default = 6.0
171   settings_lai_low_default = 1.0
172   settings_patch_height_default = 10.0
173   settings_lai_alpha = 5.0
174   settings_lai_beta = 3.0
175   ndomains = 0
176
177   global_acronym = ""
178   global_angle = ""
179   global_author = ""
180   global_campaign = ""
181   global_comment = ""
182   global_contact = ""
183   global_data_content = ""
184   global_dependencies = ""
185   global_institution = ""
186   global_keywords = ""
187   global_location = ""
188   global_palm_version = 6.0
189   global_references = ""
190   global_site = ""
191   global_source = ""
192
193   path_out = ""
194   version_out = 1
195   filename_out = "default"
196   
197   domain_names = []
198   domain_px = []
199   domain_x0 = []
200   domain_y0 = []
201   domain_x1 = []
202   domain_y1 = []
203   domain_nx = []
204   domain_ny = []
205   domain_dz = []
206   domain_3d = []
207   domain_high_vegetation = []
208   domain_ip = []
209   domain_za = []
210   domain_parent = []
211   domain_green_roofs = []
212   domain_street_trees = []
213   domain_canopy_patches = []
214
215   zt_min = 0.0
216   zt_all = []
217
218   input_names = []
219   input_px = []
220
221   input_file_x = []
222   input_file_y = [] 
223   input_file_x_UTM = []
224   input_file_y_UTM = [] 
225   input_file_lat = []
226   input_file_lon = [] 
227   
228   input_file_zt = []
229   input_file_buildings_2d = []
230   input_file_bridges_2d = []
231   input_file_building_id = []
232   input_file_bridges_id = []   
233   input_file_building_type = []
234   input_file_lai = []
235   input_file_vegetation_type = []
236   input_file_vegetation_height = []   
237   input_file_pavement_type = []
238   input_file_water_type = []
239   input_file_street_type = [] 
240   input_file_street_crossings = []     
241   input_file_soil_type = []
242   input_file_vegetation_on_roofs = []
243   input_file_tree_crown_diameter = []
244   input_file_tree_height = []
245   input_file_tree_trunk_diameter = []
246   input_file_tree_type = []
247   input_file_patch_height = []
248   
249# Load all user parameters from config file
250   for i in range(0,len(config.sections())):
251
252      read_tmp = config.sections()[i]
253     
254      if ( read_tmp == 'global' ):
255         global_acronym = config.get(read_tmp, 'acronym')
256         global_angle = float(config.get(read_tmp, 'rotation_angle'))
257         global_author = config.get(read_tmp, 'author')
258         global_campaign = config.get(read_tmp, 'campaign')
259         global_comment = config.get(read_tmp, 'comment')
260         global_contact = config.get(read_tmp, 'contact_person')
261         global_data_content = config.get(read_tmp, 'data_content')   
262         global_dependencies = config.get(read_tmp, 'dependencies')
263         global_institution = config.get(read_tmp, 'institution')
264         global_keywords = config.get(read_tmp, 'keywords')
265         global_location = config.get(read_tmp, 'location')
266         global_palm_version = float(config.get(read_tmp, 'palm_version'))
267         global_references = config.get(read_tmp, 'references')
268         global_site = config.get(read_tmp, 'site')
269         global_source = config.get(read_tmp, 'source')
270
271
272      if ( read_tmp == 'settings' ):
273         settings_lai_roof_intensive = config.get(read_tmp, 'lai_roof_intensive')
274         settings_lai_roof_extensive = config.get(read_tmp, 'lai_roof_extensive') 
275         settings_bridge_width = float(config.get(read_tmp, 'bridge_width'))   
276         settings_season = config.get(read_tmp, 'season')   
277         settings_lai_high_default = float(config.get(read_tmp, 'lai_high_vegetation_default'))
278         settings_lai_low_default = float(config.get(read_tmp, 'lai_low_vegetation_default'))         
279         settings_patch_height_default = float(config.get(read_tmp, 'patch_height_default')) 
280         settings_lai_alpha = float(config.get(read_tmp, 'lai_alpha')) 
281         settings_lai_beta = float(config.get(read_tmp, 'lai_beta')) 
282
283      if ( read_tmp == 'output' ):
284         path_out = config.get(read_tmp, 'path')         
285         filename_out = config.get(read_tmp, 'file_out')
286         version_out = float(config.get(read_tmp, 'version'))
287
288      if ( read_tmp.split("_")[0] == 'domain' ):
289         ndomains = ndomains + 1
290         domain_names.append(read_tmp.split("_")[1])
291         domain_px.append(float(config.get(read_tmp, 'pixel_size')))
292         domain_nx.append(int(config.get(read_tmp, 'nx')))
293         domain_ny.append(int(config.get(read_tmp, 'ny')))
294         domain_dz.append(float(config.get(read_tmp, 'dz')))       
295         domain_3d.append(config.getboolean(read_tmp, 'buildings_3d'))
296         domain_high_vegetation.append(config.getboolean(read_tmp, 'allow_high_vegetation'))
297         domain_canopy_patches.append(config.getboolean(read_tmp, 'generate_vegetation_patches'))       
298         domain_ip.append(config.getboolean(read_tmp, 'interpolate_terrain'))
299         domain_za.append(config.getboolean(read_tmp, 'use_palm_z_axis')) 
300         if domain_ip[ndomains-1] and not domain_za[ndomains-1]:
301            domain_za[ndomains-1] = True
302            print("+++ Overwrite user setting for use_palm_z_axis")
303         
304         domain_parent.append(config.get(read_tmp, 'domain_parent'))
305
306         domain_x0.append(int(floor(float(config.get(read_tmp, 'origin_x'))/float(config.get(read_tmp, 'pixel_size')))))
307         domain_y0.append(int(floor(float(config.get(read_tmp, 'origin_y'))/float(config.get(read_tmp, 'pixel_size')))))       
308         domain_x1.append(int(floor(float(config.get(read_tmp, 'origin_x'))/float(config.get(read_tmp, 'pixel_size'))) + int(config.get(read_tmp, 'nx'))))
309         domain_y1.append(int(floor(float(config.get(read_tmp, 'origin_y'))/float(config.get(read_tmp, 'pixel_size'))) + int(config.get(read_tmp, 'ny'))))
310         domain_green_roofs.append(config.getboolean(read_tmp, 'vegetation_on_roofs')) 
311         domain_street_trees.append(config.getboolean(read_tmp, 'street_trees')) 
312         
313      if ( read_tmp.split("_")[0] == 'input' ):
314         input_names.append(read_tmp.split("_")[1])
315         input_px.append(float(config.get(read_tmp, 'pixel_size')))
316         input_file_x.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_x'))
317         input_file_y.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_y')) 
318         input_file_lat.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_lat'))
319         input_file_lon.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_lon'))   
320         input_file_x_UTM.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_x_UTM'))
321         input_file_y_UTM.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_y_UTM')) 
322         input_file_zt.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_zt'))
323         input_file_buildings_2d.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_buildings_2d'))
324         input_file_bridges_2d.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_bridges_2d'))
325         input_file_building_id.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_building_id')) 
326         input_file_bridges_id.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_bridges_id')) 
327         input_file_building_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_building_type'))
328         input_file_lai.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_lai'))         
329         input_file_vegetation_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_vegetation_type')) 
330         input_file_vegetation_height.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_vegetation_height'))         
331         input_file_pavement_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_pavement_type'))   
332         input_file_water_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_water_type'))
333         input_file_street_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_street_type'))   
334         input_file_street_crossings.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_street_crossings')) 
335         input_file_patch_height.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_patch_height')) 
336         
337         tmp = config.get(read_tmp, 'file_tree_crown_diameter')
338         if tmp is not None:
339            input_file_tree_crown_diameter.append(config.get(read_tmp, 'path') + "/" + tmp) 
340         else:
341            input_file_tree_crown_diameter.append(None)
342         input_file_tree_height.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_tree_height')) 
343         input_file_tree_trunk_diameter.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_tree_trunk_diameter')) 
344         input_file_tree_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_tree_type')) 
345         input_file_vegetation_on_roofs.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_vegetation_on_roofs'))           
346         #input_file_soil_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_soil_type'))   
347   return 0
348
349
350############################################################
351
352# Start of main program
353datatypes = {
354   "x": "f4", 
355   "y": "f4", 
356   "z": "f4", 
357   "lat": "f4", 
358   "lon": "f4",   
359   "E_UTM": "f4", 
360   "N_UTM": "f4", 
361   "zt": "f4",
362   "buildings_2d": "f4",
363   "buildings_3d": "b",
364   "bridges_2d": "f4",
365   "building_id": "i",
366   "bridges_id": "i",
367   "building_type": "b",
368   "nsurface_fraction": "i",
369   "vegetation_type": "b",
370   "vegetation_height": "f4",
371   "pavement_type": "b", 
372   "water_type": "b", 
373   "street_type": "b", 
374   "street_crossings": "b",     
375   "soil_type": "b",
376   "surface_fraction": "f4",
377   "building_pars": "f4",
378   "vegetation_pars": "f4",
379   "tree_data": "f4",
380   "tree_type": "b",
381   "nbuilding_pars": "i",
382   "nvegetation_pars": "i",
383   "zlad": "f4"
384   }
385
386fillvalues = {
387   "lat": float(-9999.0),
388   "lon": float(-9999.0),
389   "E_UTM": float(-9999.0),
390   "N_UTM": float(-9999.0),   
391   "zt": float(-9999.0),
392   "buildings_2d": float(-9999.0),
393   "buildings_3d": np.byte(-127),
394   "bridges_2d": float(-9999.0),
395   "building_id": int(-9999),
396   "bridges_id": int(-9999),
397   "building_type": np.byte(-127),
398   "nsurface_fraction": int(-9999),
399   "vegetation_type": np.byte(-127),
400   "vegetation_height": float(-9999.0),
401   "pavement_type": np.byte(-127),
402   "water_type": np.byte(-127),
403   "street_type": np.byte(-127), 
404   "street_crossings": np.byte(-127),   
405   "soil_type": np.byte(-127),
406   "surface_fraction": float(-9999.0),
407   "building_pars": float(-9999.0),
408   "vegetation_pars": float(-9999.0),
409   "tree_data": float(-9999.0),
410   "tree_type": np.byte(-127)
411   }
412
413defaultvalues = {
414   "lat": float(-9999.0),
415   "lon": float(-9999.0),
416   "E_UTM": float(-9999.0),
417   "N_UTM": float(-9999.0), 
418   "zt": float(0.0),
419   "buildings_2d": float(0.0),
420   "buildings_3d": 0,   
421   "bridges_2d": float(0.0),
422   "building_id": int(0),
423   "bridges_id": int(0),
424   "building_type": 1,
425   "nsurface_fraction": int(-9999),
426   "vegetation_type": 3,
427   "vegetation_height": float(-9999.0),
428   "pavement_type": 1, 
429   "water_type": 1, 
430   "street_type": 1, 
431   "street_crossings": 0,   
432   "soil_type": 1,
433   "surface_fraction": float(0.0),
434   "buildings_pars": float(-9999.0),
435   "tree_data": float(-9999.0),
436   "tree_type": 0,
437   "vegetation_pars": float(-9999.0)   
438   }
439 
440# Read configuration file and set parameters accordingly
441read_config_file()
442
443
444filename = []
445ii = []
446ii_parent = []
447# Define indices and filenames for all domains and create netCDF files
448for i in range(0,ndomains):
449
450#  Calculate indices and input files
451   ii.append(input_px.index(domain_px[i]))
452   filename.append(path_out + "/" + filename_out + "_" + domain_names[i])
453   if domain_parent[i] is not None:
454      ii_parent.append(domain_names.index(domain_parent[i]))
455   else:
456      ii_parent.append(None)
457
458
459   x_UTM = nc_read_from_file_2d(input_file_x[ii[i]], "Band1", domain_x0[i], domain_x0[i], domain_y0[i], domain_y0[i])
460   y_UTM = nc_read_from_file_2d(input_file_y[ii[i]], "Band1", domain_x0[i], domain_x0[i], domain_y0[i], domain_y0[i])
461   lat = nc_read_from_file_2d(input_file_lat[ii[i]], "Band1", domain_x0[i], domain_x0[i], domain_y0[i], domain_y0[i])
462   lon = nc_read_from_file_2d(input_file_lon[ii[i]], "Band1", domain_x0[i], domain_x0[i], domain_y0[i], domain_y0[i]) 
463   
464#  Create NetCDF output file and set global attributes
465   nc_create_file(filename[i])
466   nc_write_global_attributes(filename[i],float(x_UTM[0,0]),float(y_UTM[0,0]),float(lat[0,0]),float(lon[0,0]),"",global_acronym,global_angle,global_author,global_campaign,global_comment,global_contact,global_data_content,global_dependencies,global_institution,global_keywords,global_location,global_palm_version,global_references,global_site,global_source,version_out)
467
468   
469# Process terrain height
470for i in range(0,ndomains):
471
472#  Read and write terrain height (zt)
473   zt = nc_read_from_file_2d(input_file_zt[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])
474
475#  Final step: add zt array to the global array   
476   zt_all.append(zt)
477   del zt
478 
479#  Calculate the global (all domains) minimum of the terrain height. This value will be substracted for all
480#  data sets
481zt_min = min(zt_all[0].flatten()) 
482for i in range(0,ndomains):
483   zt_min = min(zt_min,min(zt_all[i].flatten()))
484   
485del zt_all[:]
486
487print( "Shift terrain height by -" + str(zt_min))
488for i in range(0,ndomains):
489
490#  Read and write terrain height (zt)
491   zt = nc_read_from_file_2d(input_file_zt[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])
492   x = nc_read_from_file_1d(input_file_x[ii[i]], "x", domain_x0[i], domain_x1[i])
493   y = nc_read_from_file_1d(input_file_y[ii[i]], "y", domain_y0[i], domain_y1[i])
494
495   zt = zt - zt_min
496   
497   nc_write_global_attribute(filename[i], 'origin_z', float(zt_min))
498
499#  If necessary, interpolate parent domain terrain height on child domain grid and blend the two
500   if domain_ip[i]:
501      tmp_x0 = int( domain_x0[i] * domain_px[i] / domain_px[ii_parent[i]] ) - 1
502      tmp_y0 = int( domain_y0[i] * domain_px[i] / domain_px[ii_parent[i]] ) - 1
503      tmp_x1 = int( domain_x1[i] * domain_px[i] / domain_px[ii_parent[i]] ) + 1
504      tmp_y1 = int( domain_y1[i] * domain_px[i] / domain_px[ii_parent[i]] ) + 1
505               
506      tmp_x = nc_read_from_file_1d(input_file_x[ii_parent[i]], "x", tmp_x0, tmp_x1)
507      tmp_y = nc_read_from_file_1d(input_file_y[ii_parent[i]], "y", tmp_y0, tmp_y1)
508
509      zt_parent = nc_read_from_file_2d(input_file_zt[ii_parent[i]], 'Band1', tmp_x0, tmp_x1, tmp_y0, tmp_y1)
510
511      print( "Shift terrain height by -" + str(zt_min))
512      zt_parent = zt_parent - zt_min
513
514#     Interpolate array and bring to PALM grid of child domain
515
516      zt_ip = interpolate_2d(zt_parent,tmp_x,tmp_y,x,y)
517      zt_ip = bring_to_palm_grid(zt_ip,x,y,domain_dz[ii_parent[i]])
518
519#     Shift the child terrain height according to the parent mean terrain height       
520      zt = zt - np.mean(zt)  + np.mean(zt_ip)
521 
522 
523#     Blend over the parent and child terrain height within a radius of 50 px 
524      zt = blend_array_2d(zt,zt_ip,50)
525   
526#  Final step: add zt array to the global array   
527   zt_all.append(zt)
528   del zt
529 
530
531# Read and shift x and y coordinates, shift terrain height according to its minimum value and write all data
532# to file 
533for i in range(0,ndomains):   
534#  Read horizontal grid variables from zt file and write them to output file
535   x = nc_read_from_file_1d(input_file_x[ii[i]], "x", domain_x0[i], domain_x1[i])
536   y = nc_read_from_file_1d(input_file_y[ii[i]], "y", domain_y0[i], domain_y1[i]) 
537   x = x - min(x.flatten()) + domain_px[i]/2.0
538   y = y - min(y.flatten()) + domain_px[i]/2.0 
539   nc_write_dimension(filename[i], 'x', x, datatypes["x"])
540   nc_write_dimension(filename[i], 'y', y, datatypes["y"]) 
541   nc_write_attribute(filename[i], 'x', 'long_name', 'x')
542   nc_write_attribute(filename[i], 'x', 'standard_name','projection_x_coordinate')   
543   nc_write_attribute(filename[i], 'x', 'units', 'm')
544   nc_write_attribute(filename[i], 'y', 'long_name', 'x')
545   nc_write_attribute(filename[i], 'y', 'standard_name', 'projection_y_coordinate')   
546   nc_write_attribute(filename[i], 'y', 'units', 'm')
547
548   lat = nc_read_from_file_2d(input_file_lat[ii[i]], "Band1", domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])
549   lon = nc_read_from_file_2d(input_file_lon[ii[i]], "Band1", domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) 
550
551   nc_write_to_file_2d(filename[i], 'lat', lat, datatypes["lat"],'y','x',fillvalues["lat"])
552   nc_write_to_file_2d(filename[i], 'lon', lon, datatypes["lon"],'y','x',fillvalues["lon"]) 
553   
554   nc_write_attribute(filename[i], 'lat', 'long_name', 'latitude')
555   nc_write_attribute(filename[i], 'lat', 'standard_name','latitude')   
556   nc_write_attribute(filename[i], 'lat', 'units', 'degrees_north')
557 
558   nc_write_attribute(filename[i], 'lon', 'long_name', 'longitude')
559   nc_write_attribute(filename[i], 'lon', 'standard_name','longitude')   
560   nc_write_attribute(filename[i], 'lon', 'units', 'degrees_east') 
561   
562   x_UTM = nc_read_from_file_2d(input_file_x_UTM[ii[i]], "Band1", domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])
563   y_UTM = nc_read_from_file_2d(input_file_y_UTM[ii[i]], "Band1", domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) 
564   
565   nc_write_to_file_2d(filename[i], 'E_UTM', x_UTM, datatypes["E_UTM"],'y','x',fillvalues["E_UTM"])
566   nc_write_to_file_2d(filename[i], 'N_UTM', y_UTM, datatypes["N_UTM"],'y','x',fillvalues["N_UTM"]) 
567   
568   nc_write_attribute(filename[i], 'E_UTM', 'long_name', 'easting')
569   nc_write_attribute(filename[i], 'E_UTM', 'standard_name','projection_x_coorindate')   
570   nc_write_attribute(filename[i], 'E_UTM', 'units', 'm')
571 
572   nc_write_attribute(filename[i], 'N_UTM', 'long_name', 'northing')
573   nc_write_attribute(filename[i], 'N_UTM', 'standard_name','projection_y_coorindate')   
574   nc_write_attribute(filename[i], 'N_UTM', 'units', 'm') 
575 
576   nc_write_crs(filename[i])
577   
578 
579 
580#  If necessary, bring terrain height to PALM's vertical grid. This is either forced by the user or implicitly
581#  by using interpolation for a child domain       
582   if domain_za[i]:
583      zt_all[i] = bring_to_palm_grid(zt_all[i],x,y,domain_dz[i])
584
585   nc_write_to_file_2d(filename[i], 'zt', zt_all[i], datatypes["zt"],'y','x',fillvalues["zt"])
586   nc_write_attribute(filename[i], 'zt', 'long_name', 'orography')
587   nc_write_attribute(filename[i], 'zt', 'units', 'm')
588   nc_write_attribute(filename[i], 'zt', 'res_orig', domain_px[i]) 
589   nc_write_attribute(filename[i], 'zt', 'coordinates', 'E_UTM N_UTM lon lat') 
590   nc_write_attribute(filename[i], 'zt', 'grid_mapping', 'E_UTM N_UTM lon lat')   
591
592del zt_all
593
594
595#  Process building height, id, and type
596for i in range(0,ndomains): 
597   buildings_2d = nc_read_from_file_2d(input_file_buildings_2d[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])
598   
599   building_id = nc_read_from_file_2d(input_file_building_id[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) 
600 
601   building_type = nc_read_from_file_2d(input_file_building_type[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) 
602   building_type[building_type >= 254] = fillvalues["building_type"]
603   building_type = np.where(building_type < 1,defaultvalues["building_type"],building_type)
604 
605   check = check_arrays_2(buildings_2d,building_id,fillvalues["buildings_2d"],fillvalues["building_id"])
606   if not check:
607      buildings_2d = np.where(building_id != fillvalues["building_id"],buildings_2d,fillvalues["buildings_2d"])
608      building_id = np.where(buildings_2d == fillvalues["buildings_2d"],fillvalues["building_id"],building_id)
609      print("Data check #1 " + str(check_arrays_2(buildings_2d,building_id,fillvalues["buildings_2d"],fillvalues["building_id"])))
610
611   check = check_arrays_2(buildings_2d,building_type,fillvalues["buildings_2d"],fillvalues["building_type"])
612   if not check:
613      building_type = np.where(buildings_2d == fillvalues["buildings_2d"],fillvalues["building_type"],building_type)
614      building_type = np.where((building_type == fillvalues["building_type"]) & (buildings_2d != fillvalues["buildings_2d"]),defaultvalues["building_type"],building_type)
615      print("Data check #2 " + str(check_arrays_2(buildings_2d,building_type,fillvalues["buildings_2d"],fillvalues["building_type"])))
616 
617 
618   nc_write_to_file_2d(filename[i], 'buildings_2d', buildings_2d, datatypes["buildings_2d"],'y','x',fillvalues["buildings_2d"])
619   nc_write_attribute(filename[i], 'buildings_2d', 'long_name', 'buildings')
620   nc_write_attribute(filename[i], 'buildings_2d', 'units', 'm')
621   nc_write_attribute(filename[i], 'buildings_2d', 'res_orig', domain_px[i]) 
622   nc_write_attribute(filename[i], 'buildings_2d', 'lod', 1) 
623   nc_write_attribute(filename[i], 'buildings_2d', 'coordinates', 'E_UTM N_UTM lon lat') 
624   nc_write_attribute(filename[i], 'buildings_2d', 'grid_mapping', 'E_UTM N_UTM lon lat')   
625   
626   nc_write_to_file_2d(filename[i], 'building_id', building_id, datatypes["building_id"],'y','x',fillvalues["building_id"])
627   nc_write_attribute(filename[i], 'building_id', 'long_name', 'building id')
628   nc_write_attribute(filename[i], 'building_id', 'units', '')
629   nc_write_attribute(filename[i], 'building_id', 'res _orig', domain_px[i]) 
630   nc_write_attribute(filename[i], 'building_id', 'coordinates', 'E_UTM N_UTM lon lat') 
631   nc_write_attribute(filename[i], 'building_id', 'grid_mapping', 'E_UTM N_UTM lon lat')   
632   
633   nc_write_to_file_2d(filename[i], 'building_type', building_type, datatypes["building_type"],'y','x',fillvalues["building_type"])
634   nc_write_attribute(filename[i], 'building_type', 'long_name', 'building type')
635   nc_write_attribute(filename[i], 'building_type', 'units', '')
636   nc_write_attribute(filename[i], 'building_type', 'res_orig', domain_px[i]) 
637   nc_write_attribute(filename[i], 'building_type', 'coordinates', 'E_UTM N_UTM lon lat') 
638   nc_write_attribute(filename[i], 'building_type', 'grid_mapping', 'E_UTM N_UTM lon lat')   
639 
640   print("out_1: " + str(np.any(building_type == -2)))
641 
642del buildings_2d
643del building_id
644del building_type
645
646# Create 3d buildings if necessary. In that course, read bridge objects and add them to building layer
647for i in range(0,ndomains): 
648   
649   if domain_3d[i]:
650      x = nc_read_from_file_2d_all(filename[i], 'x')
651      y = nc_read_from_file_2d_all(filename[i], 'y') 
652      buildings_2d = nc_read_from_file_2d_all(filename[i], 'buildings_2d') 
653      building_id = nc_read_from_file_2d_all(filename[i], 'building_id') 
654     
655      bridges_2d = nc_read_from_file_2d(input_file_bridges_2d[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])
656      bridges_id = nc_read_from_file_2d(input_file_bridges_id[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])
657     
658      bridges_2d = np.where(bridges_2d == 0.0,fillvalues["bridges_2d"],bridges_2d)
659      building_id = np.where(bridges_2d == fillvalues["bridges_2d"],building_id,bridges_id)
660     
661     
662      if np.any(buildings_2d != fillvalues["buildings_2d"]):
663         buildings_3d, z = make_3d_from_2d(buildings_2d,x,y,domain_dz[i])
664         if np.any(bridges_2d != fillvalues["bridges_2d"]):
665            buildings_3d = make_3d_from_bridges_2d(buildings_3d,bridges_2d,x,y,domain_dz[i],settings_bridge_width,fillvalues["bridges_2d"])
666         else:
667            print("Skipping creation of 3D bridges (no bridges in domain)")
668           
669           
670         nc_write_dimension(filename[i], 'z', z, datatypes["z"]) 
671         nc_write_attribute(filename[i], 'z', 'long_name', 'z') 
672         nc_write_attribute(filename[i], 'z', 'units', 'm')
673 
674         nc_write_to_file_3d(filename[i], 'buildings_3d', buildings_3d, datatypes["buildings_3d"],'z','y','x',fillvalues["buildings_3d"])   
675         nc_write_attribute(filename[i], 'buildings_3d', 'long_name', 'buildings 3d')
676         nc_write_attribute(filename[i], 'buildings_3d', 'units', '')
677         nc_write_attribute(filename[i], 'buildings_3d', 'res_orig', domain_px[i]) 
678         nc_write_attribute(filename[i], 'buildings_3d', 'lod', 2) 
679         
680         del buildings_3d
681         
682      else:
683         print("Skipping creation of 3D buildings (no buildings in domain)")
684
685
686      del bridges_2d, bridges_id, building_id, buildings_2d
687
688
689
690# Read vegetation type, water_type, pavement_type, soil_type and make fields consistent
691for i in range(0,ndomains): 
692
693   building_type = nc_read_from_file_2d_all(filename[i], 'building_type') 
694     
695   vegetation_type = nc_read_from_file_2d(input_file_vegetation_type[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])   
696   vegetation_type[vegetation_type == 255] = fillvalues["vegetation_type"]
697   vegetation_type = np.where((vegetation_type < 1) & (vegetation_type != fillvalues["vegetation_type"]),defaultvalues["vegetation_type"],vegetation_type)
698   
699   pavement_type = nc_read_from_file_2d(input_file_pavement_type[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])   
700   pavement_type[pavement_type == 255] = fillvalues["pavement_type"]
701   pavement_type = np.where((pavement_type < 1) & (pavement_type != fillvalues["pavement_type"]),defaultvalues["pavement_type"],pavement_type)
702
703   water_type = nc_read_from_file_2d(input_file_water_type[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])   
704   water_type[water_type == 255] = fillvalues["water_type"]
705   water_type = np.where((water_type < 1) & (water_type != fillvalues["water_type"]),defaultvalues["water_type"],water_type)
706 
707#  to do: replace by real soil input data
708   soil_type = nc_read_from_file_2d(input_file_vegetation_type[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) 
709   soil_type[soil_type == 255] = fillvalues["soil_type"]
710   soil_type = np.where((soil_type < 1) & (soil_type != fillvalues["soil_type"]),defaultvalues["soil_type"],soil_type)
711
712#  Make arrays consistent
713#  #1 Set vegetation type to missing for pixel where a pavement type is set
714   vegetation_type = np.where((vegetation_type != fillvalues["vegetation_type"]) & (pavement_type != fillvalues["pavement_type"]),fillvalues["vegetation_type"],vegetation_type)
715
716#  #2 Set vegetation type to missing for pixel where a building type is set
717   vegetation_type = np.where((vegetation_type != fillvalues["vegetation_type"]) & (building_type != fillvalues["building_type"]) ,fillvalues["vegetation_type"],vegetation_type)
718
719#  #3 Set vegetation type to missing for pixel where a building type is set
720   vegetation_type = np.where((vegetation_type != fillvalues["vegetation_type"]) & (water_type != fillvalues["water_type"]),fillvalues["vegetation_type"],vegetation_type)   
721
722#  #4 Remove pavement for pixels with buildings
723   pavement_type = np.where((pavement_type != fillvalues["pavement_type"]) & (building_type != fillvalues["building_type"]),fillvalues["pavement_type"],pavement_type) 
724
725#  #5 Remove pavement for pixels with water
726   pavement_type = np.where((pavement_type != fillvalues["pavement_type"]) & (water_type != fillvalues["water_type"]),fillvalues["pavement_type"],pavement_type)   
727
728#  #6 Remove water for pixels with buildings
729   water_type = np.where((water_type != fillvalues["water_type"]) & (building_type != fillvalues["building_type"]),fillvalues["water_type"],water_type) 
730 
731   
732#  Correct vegetation_type when a vegetation height is available and is indicative of low vegeetation
733   vegetation_height = nc_read_from_file_2d(input_file_vegetation_height[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])   
734         
735   vegetation_type = np.where((vegetation_height != fillvalues["vegetation_height"]) & (vegetation_height == 0.0) & ((vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18)), 3, vegetation_type)
736   vegetation_height = np.where((vegetation_height != fillvalues["vegetation_height"]) & (vegetation_height == 0.0) & ((vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18)), fillvalues["vegetation_height"],vegetation_height)
737 
738#  Check for consistency and fill empty fields with default vegetation type
739   consistency_array, test = check_consistency_4(vegetation_type,building_type,pavement_type,water_type,fillvalues["vegetation_type"],fillvalues["building_type"],fillvalues["pavement_type"],fillvalues["water_type"])
740   
741   if test:
742      vegetation_type = np.where(consistency_array == 0,defaultvalues["vegetation_type"],vegetation_type)
743      consistency_array, test = check_consistency_4(vegetation_type,building_type,pavement_type,water_type,fillvalues["vegetation_type"],fillvalues["building_type"],fillvalues["pavement_type"],fillvalues["water_type"])
744
745#  #7 to be removed: set default soil type everywhere
746   soil_type = np.where((vegetation_type != fillvalues["vegetation_type"]) | (pavement_type != fillvalues["pavement_type"]),defaultvalues["soil_type"],fillvalues["soil_type"]) 
747       
748   
749#  Check for consistency and fill empty fields with default vegetation type
750   consistency_array, test = check_consistency_3(vegetation_type,pavement_type,soil_type,fillvalues["vegetation_type"],fillvalues["pavement_type"],fillvalues["soil_type"])
751       
752#  Create surface_fraction array
753   x = nc_read_from_file_2d_all(filename[i], 'x')
754   y = nc_read_from_file_2d_all(filename[i], 'y') 
755   nsurface_fraction = np.arange(0,3)
756   surface_fraction = np.ones((len(nsurface_fraction),len(y),len(x)))
757   
758   surface_fraction[0,:,:] = np.where(vegetation_type != fillvalues["vegetation_type"], 1.0, 0.0)
759   surface_fraction[1,:,:] = np.where(pavement_type != fillvalues["pavement_type"], 1.0, 0.0)
760   surface_fraction[2,:,:] = np.where(water_type != fillvalues["water_type"], 1.0, 0.0)
761   
762   nc_write_dimension(filename[i], 'nsurface_fraction', nsurface_fraction, datatypes["nsurface_fraction"]) 
763   nc_write_to_file_3d(filename[i], 'surface_fraction', surface_fraction, datatypes["surface_fraction"],'nsurface_fraction','y','x',fillvalues["surface_fraction"])     
764   nc_write_attribute(filename[i], 'surface_fraction', 'long_name', 'surface fraction')
765   nc_write_attribute(filename[i], 'surface_fraction', 'units', '')
766   nc_write_attribute(filename[i], 'surface_fraction', 'res_orig', domain_px[i])     
767   del surface_fraction   
768
769   nc_write_to_file_2d(filename[i], 'vegetation_type', vegetation_type, datatypes["vegetation_type"],'y','x',fillvalues["vegetation_type"])     
770   nc_write_attribute(filename[i], 'vegetation_type', 'long_name', 'vegetation type')
771   nc_write_attribute(filename[i], 'vegetation_type', 'units', '')
772   nc_write_attribute(filename[i], 'vegetation_type', 'res_orig', domain_px[i]) 
773   nc_write_attribute(filename[i], 'vegetation_type', 'coordinates', 'E_UTM N_UTM lon lat') 
774   nc_write_attribute(filename[i], 'vegetation_type', 'grid_mapping', 'E_UTM N_UTM lon lat')   
775   del vegetation_type
776
777   nc_write_to_file_2d(filename[i], 'pavement_type', pavement_type, datatypes["pavement_type"],'y','x',fillvalues["pavement_type"])   
778   nc_write_attribute(filename[i], 'pavement_type', 'long_name', 'pavement type')
779   nc_write_attribute(filename[i], 'pavement_type', 'units', '')
780   nc_write_attribute(filename[i], 'pavement_type', 'res_orig', domain_px[i]) 
781   nc_write_attribute(filename[i], 'pavement_type', 'coordinates', 'E_UTM N_UTM lon lat') 
782   nc_write_attribute(filename[i], 'pavement_type', 'grid_mapping', 'E_UTM N_UTM lon lat')   
783   del pavement_type
784
785   nc_write_to_file_2d(filename[i], 'water_type', water_type, datatypes["water_type"],'y','x',fillvalues["water_type"]) 
786   nc_write_attribute(filename[i], 'water_type', 'long_name', 'water type')
787   nc_write_attribute(filename[i], 'water_type', 'units', '')
788   nc_write_attribute(filename[i], 'water_type', 'res_orig', domain_px[i]) 
789   nc_write_attribute(filename[i], 'water_type', 'coordinates', 'E_UTM N_UTM lon lat') 
790   nc_write_attribute(filename[i], 'water_type', 'grid_mapping', 'E_UTM N_UTM lon lat')   
791   del water_type
792
793   nc_write_to_file_2d(filename[i], 'soil_type', soil_type, datatypes["soil_type"],'y','x',fillvalues["soil_type"]) 
794   nc_write_attribute(filename[i], 'soil_type', 'long_name', 'soil type')
795   nc_write_attribute(filename[i], 'soil_type', 'units', '')
796   nc_write_attribute(filename[i], 'soil_type', 'res_orig', domain_px[i])   
797   nc_write_attribute(filename[i], 'soil_type', 'coordinates', 'E_UTM N_UTM lon lat') 
798   nc_write_attribute(filename[i], 'soil_type', 'grid_mapping', 'E_UTM N_UTM lon lat')   
799   del soil_type
800
801   del x
802   del y
803
804#  pixels with bridges get building_type = 7 = bridge. This does not change the _type setting for the under-bridge
805#  area NOTE: when bridges are present the consistency check will fail at the moment
806   if domain_3d[i]:
807      if np.any(building_type != fillvalues["building_type"]):
808 
809         bridges_2d = nc_read_from_file_2d(input_file_bridges_2d[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])         
810         bridges_2d = np.where(bridges_2d == 0.0,fillvalues["bridges_2d"],bridges_2d)
811         building_type = np.where(bridges_2d != fillvalues["bridges_2d"],7,building_type)
812         nc_overwrite_to_file_2d(filename[i], 'building_type', building_type) 
813 
814         del building_type 
815         del bridges_2d
816
817# Read/write street type and street crossings
818for i in range(0,ndomains): 
819
820   street_type = nc_read_from_file_2d(input_file_street_type[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])   
821   street_type[street_type == 255] = fillvalues["street_type"]
822   street_type = np.where((street_type < 1) & (street_type != fillvalues["street_type"]),defaultvalues["street_type"],street_type)
823   
824   nc_write_to_file_2d(filename[i], 'street_type', street_type, datatypes["street_type"],'y','x',fillvalues["street_type"]) 
825   nc_write_attribute(filename[i], 'street_type', 'long_name', 'street type')
826   nc_write_attribute(filename[i], 'street_type', 'units', '')
827   nc_write_attribute(filename[i], 'street_type', 'res_orig', domain_px[i]) 
828   nc_write_attribute(filename[i], 'street_type', 'coordinates', 'E_UTM N_UTM lon lat') 
829   nc_write_attribute(filename[i], 'street_type', 'grid_mapping', 'E_UTM N_UTM lon lat')   
830   del street_type
831   
832   street_crossings = nc_read_from_file_2d(input_file_street_crossings[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])   
833   street_crossings[street_crossings == 255] = fillvalues["street_crossings"]
834   street_crossings = np.where((street_crossings < 1) & (street_crossings != fillvalues["street_crossings"]),defaultvalues["street_crossings"],street_crossings)
835   
836   nc_write_to_file_2d(filename[i], 'street_crossings', street_crossings, datatypes["street_crossings"],'y','x',fillvalues["street_crossings"]) 
837   nc_write_attribute(filename[i], 'street_crossings', 'long_name', 'street crossings')
838   nc_write_attribute(filename[i], 'street_crossings', 'units', '')
839   nc_write_attribute(filename[i], 'street_crossings', 'res_orig', domain_px[i]) 
840   nc_write_attribute(filename[i], 'street_crossings', 'coordinates', 'E_UTM N_UTM lon lat') 
841   nc_write_attribute(filename[i], 'street_crossings', 'grid_mapping', 'E_UTM N_UTM lon lat')   
842   del street_crossings
843
844
845# Read/write vegetation on roofs
846for i in range(0,ndomains): 
847   if domain_green_roofs[i]: 
848      green_roofs = nc_read_from_file_2d(input_file_vegetation_on_roofs[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])   
849      buildings_2d = nc_read_from_file_2d_all(filename[i], 'buildings_2d') 
850
851
852      x = nc_read_from_file_2d_all(filename[i], 'x')
853      y = nc_read_from_file_2d_all(filename[i], 'y') 
854      nbuilding_pars = np.arange(0,46)
855      building_pars = np.ones((len(nbuilding_pars),len(y),len(x)))
856      building_pars[:,:,:] = fillvalues["building_pars"]
857
858#     assign green fraction on roofs
859      building_pars[3,:,:] = np.where( (buildings_2d != fillvalues["buildings_2d"] ) & ( green_roofs != fillvalues["building_pars"] ),1.0,fillvalues["building_pars"])
860
861#     assign leaf area index for vegetation on roofs
862      building_pars[4,:,:] = np.where( ( buildings_2d != fillvalues["buildings_2d"] ) & ( green_roofs == 1.0 ),settings_lai_roof_intensive,fillvalues["building_pars"])
863      building_pars[4,:,:] = np.where( ( buildings_2d != fillvalues["buildings_2d"] ) & ( green_roofs == 2.0 ),settings_lai_roof_extensive,building_pars[4,:,:])
864   
865   
866      nc_write_dimension(filename[i], 'nbuilding_pars', nbuilding_pars, datatypes["nbuilding_pars"]) 
867      nc_write_to_file_3d(filename[i], 'building_pars', building_pars, datatypes["building_pars"],'nbuilding_pars','y','x',fillvalues["building_pars"]) 
868      nc_write_attribute(filename[i], 'building_pars', 'long_name', 'building_pars')
869      nc_write_attribute(filename[i], 'building_pars', 'units', '')
870      nc_write_attribute(filename[i], 'building_pars', 'res_orig', domain_px[i]) 
871      nc_write_attribute(filename[i], 'building_pars', 'coordinates', 'E_UTM N_UTM lon lat') 
872      nc_write_attribute(filename[i], 'building_pars', 'grid_mapping', 'E_UTM N_UTM lon lat')   
873         
874      del building_pars, buildings_2d, x, y
875
876
877# Read tree data and create LAD and BAD arrays using the canopy generator
878for i in range(0,ndomains): 
879   lai = nc_read_from_file_2d(input_file_lai[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])   
880     
881   vegetation_type = nc_read_from_file_2d_all(filename[i], 'vegetation_type')
882     
883   lai = np.where(vegetation_type == fillvalues["vegetation_type"],fillvalues["vegetation_pars"],lai)
884
885     
886   x = nc_read_from_file_2d_all(filename[i], 'x')
887   y = nc_read_from_file_2d_all(filename[i], 'y') 
888   nvegetation_pars = np.arange(0,12)
889   vegetation_pars = np.ones((len(nvegetation_pars),len(y),len(x)))
890   vegetation_pars[:,:,:] = fillvalues["vegetation_pars"] 
891     
892   vegetation_pars[1,:,:] = lai
893     
894
895#  Write out first version of LAI. Will later be overwritten.
896   nc_write_dimension(filename[i], 'nvegetation_pars', nvegetation_pars, datatypes["nvegetation_pars"]) 
897   nc_write_to_file_3d(filename[i], 'vegetation_pars', vegetation_pars, datatypes["vegetation_pars"],'nvegetation_pars','y','x',fillvalues["vegetation_pars"]) 
898   nc_write_attribute(filename[i], 'vegetation_pars', 'long_name', 'vegetation_pars')
899   nc_write_attribute(filename[i], 'vegetation_pars', 'units', '')
900   nc_write_attribute(filename[i], 'vegetation_pars', 'res_orig', domain_px[i]) 
901   nc_write_attribute(filename[i], 'vegetation_pars', 'coordinates', 'E_UTM N_UTM lon lat') 
902   nc_write_attribute(filename[i], 'vegetation_pars', 'grid_mapping', 'E_UTM N_UTM lon lat')   
903
904   del lai, vegetation_pars, vegetation_type
905
906# Read tree data and create LAD and BAD arrays using the canopy generator
907for i in range(0,ndomains): 
908   if domain_street_trees[i]: 
909
910      vegetation_pars =  nc_read_from_file_2d_all(filename[i], 'vegetation_pars')
911
912      lai = np.copy(vegetation_pars[1,:,:])
913
914      x = nc_read_from_file_2d_all(filename[i], 'x')
915      y = nc_read_from_file_2d_all(filename[i], 'y') 
916
917#     Save lai data as default for low and high vegetation
918      lai_low = lai
919      lai_high = lai
920   
921#     Read all tree parameters from file
922      tree_height = nc_read_from_file_2d(input_file_tree_height[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])   
923
924      if (input_file_tree_crown_diameter[ii[i]] is not None):
925         tree_crown_diameter = nc_read_from_file_2d(input_file_tree_crown_diameter[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])   
926         tree_crown_diameter = np.where( (tree_crown_diameter == 0.0) | (tree_crown_diameter == -1.0) ,fillvalues["tree_data"],tree_crown_diameter)
927      else:
928         tree_crown_diameter = np.ones((len(y),len(x)))
929         tree_crown_diameter[:,:] = fillvalues["tree_data"]
930
931
932      tree_trunk_diameter = nc_read_from_file_2d(input_file_tree_trunk_diameter[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])         
933      tree_type = nc_read_from_file_2d(input_file_tree_type[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) 
934      patch_height = nc_read_from_file_2d(input_file_patch_height[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) 
935     
936#     Remove missing values from the data. Reasonable values will be set by the tree generator   
937      tree_height = np.where( (tree_height == 0.0) | (tree_height == -1.0) ,fillvalues["tree_data"],tree_height)
938      tree_trunk_diameter = np.where( (tree_trunk_diameter == 0.0) | (tree_trunk_diameter == -1.0) ,fillvalues["tree_data"],tree_trunk_diameter)
939      tree_type = np.where( (tree_type == 0.0) | (tree_type == -1.0) ,fillvalues["tree_data"],tree_type)
940      patch_height = np.where( patch_height == -1.0 ,fillvalues["tree_data"],patch_height)     
941
942#     Convert trunk diameter from cm to m
943      tree_trunk_diameter = np.where(tree_trunk_diameter != fillvalues["tree_data"], tree_trunk_diameter * 0.01,tree_trunk_diameter)
944
945
946#     Temporarily change missing value for tree_type
947      tree_type = np.where( (tree_type == fillvalues["tree_type"]),fillvalues["tree_data"],tree_type)   
948
949#     Compare patch height array with vegetation type and correct accordingly
950      vegetation_type = nc_read_from_file_2d_all(filename[i], 'vegetation_type') 
951
952
953#     For zero-height patches, set vegetation_type to short grass and remove these pixels from the patch height array
954      vegetation_type = np.where( (patch_height == 0.0) & ( (vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18) ),3,vegetation_type)
955      patch_type = np.where( (patch_height == 0.0) & ( (vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18) ),fillvalues["tree_data"],patch_height)   
956     
957
958      max_tree_height = max(tree_height.flatten())
959      max_patch_height = max(patch_height.flatten())
960     
961      if ( (max_tree_height != fillvalues["tree_data"]) | (max_patch_height == fillvalues["tree_data"]) ):
962
963         lad, bad, tree_ids, zlad = generate_single_tree_lad(x,y,domain_dz[i],max_tree_height,max_patch_height,tree_type,tree_height,tree_crown_diameter,tree_trunk_diameter,lai,settings_season,fillvalues["tree_data"])
964 
965 
966#        Remove LAD volumes that are inside buildings
967         buildings_2d = nc_read_from_file_2d_all(filename[i], 'buildings_2d') 
968         for k in range(0,len(zlad)-1):
969 
970            lad[k,:,:] = np.where(buildings_2d == fillvalues["buildings_2d"],lad[k,:,:],fillvalues["tree_data"])
971            bad[k,:,:] = np.where(buildings_2d == fillvalues["buildings_2d"],bad[k,:,:],fillvalues["tree_data"])
972            tree_ids[k,:,:] = np.where(buildings_2d == fillvalues["buildings_2d"],tree_ids[k,:,:],fillvalues["tree_data"])     
973
974         del buildings_2d
975 
976         nc_write_dimension(filename[i], 'zlad', zlad, datatypes["tree_data"]) 
977         nc_write_to_file_3d(filename[i], 'lad', lad, datatypes["tree_data"],'zlad','y','x',fillvalues["tree_data"]) 
978         nc_write_attribute(filename[i], 'lad', 'long_name', 'leaf area density')
979         nc_write_attribute(filename[i], 'lad', 'units', '')
980         nc_write_attribute(filename[i], 'lad', 'res_orig', domain_px[i]) 
981         nc_write_attribute(filename[i], 'lad', 'coordinates', 'E_UTM N_UTM lon lat') 
982         nc_write_attribute(filename[i], 'lad', 'grid_mapping', 'E_UTM N_UTM lon lat') 
983 
984         nc_write_to_file_3d(filename[i], 'bad', bad, datatypes["tree_data"],'zlad','y','x',fillvalues["tree_data"]) 
985         nc_write_attribute(filename[i], 'bad', 'long_name', 'basal area density')
986         nc_write_attribute(filename[i], 'bad', 'units', '')
987         nc_write_attribute(filename[i], 'bad', 'res_orig', domain_px[i]) 
988         nc_write_attribute(filename[i], 'bad', 'coordinates', 'E_UTM N_UTM lon lat') 
989         nc_write_attribute(filename[i], 'bad', 'grid_mapping', 'E_UTM N_UTM lon lat')
990           
991         nc_write_to_file_3d(filename[i], 'tree_id', tree_ids, datatypes["tree_data"],'zlad','y','x',fillvalues["tree_data"]) 
992         nc_write_attribute(filename[i], 'tree_id', 'long_name', 'tree id')
993         nc_write_attribute(filename[i], 'tree_id', 'units', '')
994         nc_write_attribute(filename[i], 'tree_id', 'res_orig', domain_px[i]) 
995         nc_write_attribute(filename[i], 'tree_id', 'coordinates', 'E_UTM N_UTM lon lat') 
996         nc_write_attribute(filename[i], 'tree_id', 'grid_mapping', 'E_UTM N_UTM lon lat')           
997           
998         del lai, lad, bad, tree_ids, zlad
999
1000      del vegetation_pars, tree_height, tree_crown_diameter, tree_trunk_diameter, tree_type, patch_height, x, y
1001
1002
1003# Create vegetation patches for locations with high vegetation type
1004for i in range(0,ndomains): 
1005   if domain_canopy_patches[i]: 
1006
1007#     Load vegetation_type and lad array (at level z = 0) for re-processing
1008      vegetation_type = nc_read_from_file_2d_all(filename[i], 'vegetation_type') 
1009      lad = nc_read_from_file_3d_all(filename[i], 'lad') 
1010      zlad = nc_read_from_file_1d_all(filename[i], 'zlad') 
1011      patch_height = nc_read_from_file_2d(input_file_patch_height[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) 
1012      vegetation_pars = nc_read_from_file_3d_all(filename[i], 'vegetation_pars') 
1013      lai = vegetation_pars[1,:,:]
1014
1015
1016#     Determine all pixels that do not already have and LAD but which are high vegetation to a dummy value of 1.0 and remove all other pixels
1017      lai_high = np.where( (lad[0,:,:] == fillvalues["tree_data"]) & ( ( (vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18) ) & ( (patch_height == fillvalues["tree_data"]) | (patch_height >= domain_dz[i])) ),1.0,fillvalues["tree_data"])
1018
1019#     Now, assign either the default LAI for high vegetation or the value stored in the LAI array.
1020      lai_high = np.where( (lai_high != fillvalues["tree_data"]) & (lai == fillvalues["tree_data"]), settings_lai_high_default, lai)
1021
1022#     Define a patch height wherever it is missing, but where a high vegetation LAI was set
1023      patch_height = np.where ( (lai_high != fillvalues["tree_data"]) & (patch_height == fillvalues["tree_data"]), settings_patch_height_default, patch_height)
1024
1025#     Remove pixels where street trees were already set
1026      patch_height = np.where ( (lad[0,:,:] != fillvalues["tree_data"]), fillvalues["tree_data"], patch_height)
1027
1028#     Remove patch heights that have no lai_high value
1029      patch_height = np.where ( (lai_high == fillvalues["tree_data"]), fillvalues["tree_data"], patch_height)
1030
1031#     For missing LAI values, set either the high vegetation default or the low vegetation default   
1032      lai_high = np.where( (patch_height > 2.0) & (patch_height != fillvalues["tree_data"]) & (lai_high == fillvalues["tree_data"]),settings_lai_high_default,lai_high)
1033      lai_high = np.where( (patch_height <= 2.0) & (patch_height != fillvalues["tree_data"]) & (lai_high == fillvalues["tree_data"]),settings_lai_low_default,lai_high)
1034         
1035      if ( max(patch_height.flatten()) >= (2.0 * domain_dz[i]) ):
1036         print("    start calculating LAD (this might take some time)") 
1037         lad_patch, patch_nz, status = process_patch(domain_dz[i],patch_height,max(zlad),lai_high,settings_lai_alpha,settings_lai_beta)
1038
1039         lad[0:patch_nz+1,:,:] = np.where( (lad[0:patch_nz+1,:,:] == fillvalues["tree_data"]),lad_patch[0:patch_nz+1,:,:], lad[0:patch_nz+1,:,:] )
1040
1041#     Remove high vegetation wherever it is replaced by a leaf area density. This should effectively remove all high vegetation pixels
1042      vegetation_type = np.where((lad[0,:,:] != fillvalues["tree_data"]) & (vegetation_type != fillvalues["vegetation_type"]),3,vegetation_type)
1043             
1044#     If desired, remove all high vegetation. TODO: check if this is still necessary
1045      if domain_high_vegetation[i]:
1046         vegetation_type = np.where((vegetation_type != fillvalues["vegetation_type"]) & ( (vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18) ),3,vegetation_type)   
1047
1048 
1049#     Set default low LAI for pixels with an LAD (short grass below trees)
1050      lai_low = np.where( (lad[0,:,:] == fillvalues["tree_data"]), lai, settings_lai_low_default)
1051   
1052#     Fill low vegetation pixels without LAI set or with LAI = 0 with default value
1053      lai_low = np.where( ( (lai_low == fillvalues["tree_data"]) | (lai_low == 0.0) ) & (vegetation_type != fillvalues["vegetation_type"] ), settings_lai_low_default, lai_low)
1054
1055#     Remove lai for pixels that have no vegetation_type
1056      lai_low = np.where( vegetation_type != fillvalues["vegetation_type"], lai_low, fillvalues["tree_data"])
1057       
1058#     Overwrite lai in vegetation_parameters 
1059      vegetation_pars[1,:,:] = np.copy(lai_low)
1060      nc_overwrite_to_file_3d(filename[i], 'vegetation_pars', vegetation_pars) 
1061
1062#     Overwrite lad array
1063      nc_overwrite_to_file_3d(filename[i], 'lad', lad) 
1064           
1065      nc_overwrite_to_file_2d(filename[i], 'vegetation_type', vegetation_type) 
1066
1067     
1068      del vegetation_type, lad, lai, patch_height, vegetation_pars, zlad
1069
1070# Final consistency check
1071for i in range(0,ndomains): 
1072   vegetation_type = nc_read_from_file_2d_all(filename[i], 'vegetation_type') 
1073   pavement_type = nc_read_from_file_2d_all(filename[i], 'pavement_type') 
1074   building_type = nc_read_from_file_2d_all(filename[i], 'building_type') 
1075   water_type = nc_read_from_file_2d_all(filename[i], 'water_type')     
1076   soil_type = nc_read_from_file_2d_all(filename[i], 'soil_type') 
1077
1078#  Check for consistency and fill empty fields with default vegetation type
1079   consistency_array, test = check_consistency_4(vegetation_type,building_type,pavement_type,water_type,fillvalues["vegetation_type"],fillvalues["building_type"],fillvalues["pavement_type"],fillvalues["water_type"])
1080
1081#  Check for consistency and fill empty fields with default vegetation type
1082   consistency_array, test = check_consistency_3(vegetation_type,pavement_type,soil_type,fillvalues["vegetation_type"],fillvalues["pavement_type"],fillvalues["soil_type"])
1083
1084   surface_fraction = nc_read_from_file_3d_all(filename[i], 'surface_fraction')       
1085   surface_fraction[0,:,:] = np.where(vegetation_type != fillvalues["vegetation_type"], 1.0, 0.0)
1086   surface_fraction[1,:,:] = np.where(pavement_type != fillvalues["pavement_type"], 1.0, 0.0)
1087   surface_fraction[2,:,:] = np.where(water_type != fillvalues["water_type"], 1.0, 0.0)
1088   nc_overwrite_to_file_3d(filename[i], 'surface_fraction', surface_fraction) 
1089       
1090   del vegetation_type, pavement_type, building_type, water_type, soil_type
Note: See TracBrowser for help on using the repository browser.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |