Home


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

source: palm/trunk/SCRIPTS/palm_csd @ 3756

Last change on this file since 3756 was 3726, checked in by maronga, 6 years ago

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