Index: /palm/trunk/SCRIPTS/.csd.config.default
===================================================================
--- /palm/trunk/SCRIPTS/.csd.config.default	(revision 3567)
+++ /palm/trunk/SCRIPTS/.csd.config.default	(revision 3567)
@@ -0,0 +1,171 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+#--------------------------------------------------------------------------------#
+# This file is part of the PALM model system.
+#
+# PALM is free software: you can redistribute it and/or modify it under the terms
+# of the GNU General Public License as published by the Free Software Foundation,
+# either version 3 of the License, or (at your option) any later version.
+#
+# PALM is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+# A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along with
+# PALM. If not, see <http://www.gnu.org/licenses/>.
+#
+# Copyright 1997-2018  Leibniz Universitaet Hannover
+#--------------------------------------------------------------------------------#
+#
+# Current revisions:
+# -----------------
+# 
+# 
+# Former revisions:
+# -----------------
+# $Id$
+# Initial revisions
+#
+#
+# 
+#
+#
+# Description:
+# ------------
+# Default configuration file for palm_csd
+#
+# @Author Bjoern Maronga (maronga@muk.uni-hannover.de)
+#------------------------------------------------------------------------------#
+
+[global]
+author = Bjoern Maronga, maronga@muk.uni-hannover.de
+contact_person = Bjoern Maronga, maronga@muk.uni-hannover.de
+acronym = LUHimuk
+comment = created with palm_csd
+data_content = text
+dependencies = text
+keywords = text
+source = text
+campaign = text
+location = Berlin
+site = Berlin Mitte
+institution = Institute of Meteorology and Climatology, Leibniz University Hannover
+palm_version = 6.0
+version = 1
+rotation_angle = 0.0
+references = text
+
+
+[settings]
+filename_out = /ldata2/MOSAIK/showcase_python
+lai_season = spring
+lai_high_vegetation_default = 5.0
+lai_low_vegetation_default = 1.0
+lai_alpha = 5.0
+lai_beta = 3.0
+tree_height_default = 10.0
+bridge_width = 3.0
+debug_mode = False
+
+
+[input_01]
+path = /ldata2/MOSAIK/Berlin_static_driver_data
+pixel_size = 15.0
+file_x = Berlin_CoordinatesUTM_x_15m_DLR.nc
+file_y = Berlin_CoordinatesUTM_y_15m_DLR.nc
+file_x_UTM = Berlin_CoordinatesUTM_y_15m_DLR.nc
+file_y_UTM = Berlin_CoordinatesUTM_x_15m_DLR.nc
+file_lat = Berlin_CoordinatesLatLon_y_15m_DLR.nc
+file_lon = Berlin_CoordinatesLatLon_x_15m_DLR.nc
+file_zt = Berlin_terrain_height_15m_DLR.nc
+file_buildings_2d = Berlin_building_height_15m_LOD2_DLR.nc
+file_building_id = Berlin_building_id_15m_DLR.nc
+file_building_type = Berlin_building_type_15m_DLR.nc
+file_bridges_2d = Berlin_bridges_height_15m_DLR.nc
+file_bridges_id = Berlin_bridges_id_15m_DLR.nc
+file_vegetation_type = Berlin_vegetation_type_15m_DLR.nc
+file_vegetation_height = Berlin_vegetation_patch_height_15m_DLR.nc
+file_pavement_type = Berlin_pavement_type_15m_DLR.nc
+file_water_type = Berlin_water_type_15m_DLR.nc
+file_soil_type = Berlin_soil_type_15m_DLR.nc
+file_street_type =  Berlin_street_type_15m_DLR.nc
+file_street_crossings =  Berlin_street_crossings_15m_DLR.nc
+
+[input_02]
+path = /ldata2/MOSAIK/Berlin_static_driver_data
+pixel_size = 10.0
+file_x = Berlin_CoordinatesUTM_x_10m_DLR.nc
+file_y = Berlin_CoordinatesUTM_y_10m_DLR.nc
+file_x_UTM = Berlin_CoordinatesUTM_y_10m_DLR.nc
+file_y_UTM = Berlin_CoordinatesUTM_x_10m_DLR.nc
+file_lat = Berlin_CoordinatesLatLon_y_10m_DLR.nc
+file_lon = Berlin_CoordinatesLatLon_x_10m_DLR.nc
+file_zt = Berlin_terrain_height_10m_DLR.nc
+file_buildings_2d = Berlin_building_height_10m_LOD2_DLR.nc
+file_building_id = Berlin_building_id_10m_DLR.nc
+file_building_type = Berlin_building_type_10m_DLR.nc
+file_bridges_2d = Berlin_bridges_height_10m_DLR.nc
+file_bridges_id = Berlin_bridges_id_10m_DLR.nc
+file_vegetation_type = Berlin_vegetation_type_10m_DLR.nc
+file_vegetation_height = Berlin_vegetation_patch_height_10m_DLR.nc
+file_pavement_type = Berlin_pavement_type_10m_DLR.nc
+file_water_type = Berlin_water_type_10m_DLR.nc
+file_soil_type = Berlin_soil_type_10m_DLR.nc
+file_street_type =  Berlin_street_type_10m_DLR.nc
+file_street_crossings =  Berlin_street_crossings_10m_DLR.nc
+
+[input_03]
+file_x = Berlin_CoordinatesUTM_x_1m_DLR.nc
+file_y = Berlin_CoordinatesUTM_y_1m_DLR.nc
+file_x_UTM = Berlin_CoordinatesUTM_x_1m_DLR.nc
+file_y_UTM = Berlin_CoordinatesUTM_y_1m_DLR.nc
+file_lat = Berlin_CoordinatesLatLon_y_1m_DLR.nc
+file_lon = Berlin_CoordinatesLatLon_x_1m_DLR.nc
+path = /ldata2/MOSAIK/Berlin_static_driver_data
+pixel_size = 1.0
+file_zt = Berlin_terrain_height_1m_DLR.nc
+file_buildings_2d = Berlin_building_height_1m_LOD2_DLR.nc
+file_building_id = Berlin_building_id_1m_DLR.nc
+file_bridges_2d = Berlin_bridges_height_1m_DLR.nc
+file_bridges_id = Berlin_bridges_id_1m_DLR.nc
+file_building_type = Berlin_building_type_1m_DLR.nc
+file_vegetation_type = Berlin_vegetation_type_1m_DLR.nc
+file_vegetation_height = Berlin_vegetation_patch_height_1m_DLR.nc
+file_pavement_type = Berlin_pavement_type_1m_DLR.nc
+file_water_type = Berlin_water_type_1m_DLR.nc
+file_soil_type = Berlin_soil_type_1m_DLR.nc
+file_street_type =  Berlin_street_type_1m_DLR.nc
+file_street_crossings =  Berlin_street_crossings_1m_DLR.nc
+
+[output]
+path = /ldata2/MOSAIK/
+file_out = winter_iop1_test
+version = 1
+
+[domain_root]
+pixel_size = 15.0
+origin_x = 0.0
+origin_y = 0.0
+nx = 311
+ny = 257
+buildings_3d = False
+dz = 15.0
+allow_high_vegetation = True
+generate_vegetation_patches = True
+use_palm_z_axis= True
+interpolate_terrain = False
+domain_parent
+
+[domain_N02]
+pixel_size = 1.0
+origin_x = 19605
+origin_y = 20895
+nx = 1023
+ny = 1023
+buildings_3d = True
+dz = 1.0
+allow_high_vegetation = False
+generate_vegetation_patches = True
+use_palm_z_axis= True
+interpolate_terrain = True
+domain_parent = root
Index: /palm/trunk/SCRIPTS/palm_csd
===================================================================
--- /palm/trunk/SCRIPTS/palm_csd	(revision 3567)
+++ /palm/trunk/SCRIPTS/palm_csd	(revision 3567)
@@ -0,0 +1,764 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+#--------------------------------------------------------------------------------#
+# This file is part of the PALM model system.
+#
+# PALM is free software: you can redistribute it and/or modify it under the terms
+# of the GNU General Public License as published by the Free Software Foundation,
+# either version 3 of the License, or (at your option) any later version.
+#
+# PALM is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+# A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along with
+# PALM. If not, see <http://www.gnu.org/licenses/>.
+#
+# Copyright 1997-2018  Leibniz Universitaet Hannover
+#--------------------------------------------------------------------------------#
+#
+# Current revisions:
+# -----------------
+# 
+# 
+# Former revisions:
+# -----------------
+# $Id$
+# Initial revisions
+#
+#
+# 
+#
+#
+# Description:
+# ------------
+# Processing tool for creating PIDS conform static drivers from rastered NetCDF 
+# input
+#
+# @Author Bjoern Maronga (maronga@muk.uni-hannover.de)
+#
+# @todo Remove high vegetation on demand
+# @todo Add vegetation_pars (LAI)
+# @todo Add building_pars (green roofs)
+# @todo Add LAD and BAD arrays (canopy generator)
+# @todo Make input files optional
+# @todo Allow for ASCII input of terrain height and building height
+# @todo Modularize reading config file
+#------------------------------------------------------------------------------#
+
+from palm_csd_files.palm_csd_netcdf_interface import *
+from palm_csd_files.palm_csd_tools import *
+import numpy as np
+
+    
+def read_config_file():
+   
+   import configparser
+   from math import floor
+   import numpy as np
+   import os
+   import subprocess as sub
+   import sys
+
+#  Check if configuration files exists and quit otherwise
+   input_config = ".csd.config"
+   for i in range(1,len(sys.argv)): 
+      input_config = str(sys.argv[i])
+        
+   config = configparser.RawConfigParser(allow_no_value=True)
+
+   if ( os.path.isfile(input_config) == False ):
+      print ("Error. No configuration file " + input_config + " found.")
+      raise SystemExit   
+   else:
+      print(os.path.isfile(input_config))
+
+   config.read(input_config)
+
+
+#  Definition of settings
+   global settings_filename_out
+   global settings_lai_season
+   global settings_bridge_width
+   global ndomains
+
+#  Definition of global configuration parameters
+   global global_acronym
+   global global_angle
+   global global_author
+   global global_campaign
+   global global_comment
+   global global_contact
+   global global_data_content  
+   global global_dependencies
+   global global_institution
+   global global_keywords
+   global global_location
+   global global_palm_version
+   global global_references
+   global global_site
+   global global_source
+   global global_version
+
+
+#  Definition of domain parameters
+   global domain_names
+   global domain_px
+   global domain_x0
+   global domain_y0
+   global domain_x1
+   global domain_y1
+   global domain_nx
+   global domain_ny
+   global domain_dz
+   global domain_3d
+   global domain_hv
+   global domain_cg
+   global domain_ip
+   global domain_za
+   global domain_parent
+
+#  Definition of input data parameters
+   global input_names
+   global input_px
+  
+        
+   global input_file_x
+   global input_file_y 
+   global input_file_x_UTM
+   global input_file_y_UTM   
+   global input_file_lat
+   global input_file_lon  
+   global input_file_zt
+   global input_file_buildings_2d
+   global input_file_bridges_2d 
+   global input_file_building_id
+   global input_file_bridges_id   
+   global input_file_building_type
+   global input_file_building_type
+   global input_file_vegetation_type
+   global input_file_vegetation_height
+   global input_file_pavement_type
+   global input_file_water_type
+   global input_file_street_type
+   global input_file_street_crossings
+   global input_file_soil_type
+
+
+   global zt_all
+   global zt_min
+
+   settings_filename_out = "default"
+   settings_lai_season = "summer"
+   settings_bridge_width = 3.0
+   ndomains = 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 = 6.0
+   global_references = ""
+   global_site = ""
+   global_source = ""
+   global_version = 1
+
+   domain_names = []
+   domain_px = []
+   domain_x0 = []
+   domain_y0 = []
+   domain_x1 = []
+   domain_y1 = []
+   domain_nx = []
+   domain_ny = []
+   domain_dz = []
+   domain_3d = []
+   domain_hv = []
+   domain_cg = []
+   domain_ip = []
+   domain_za = []
+   domain_parent = []
+
+   zt_min = 0.0
+   zt_all = []
+
+   input_names = []
+   input_px = []
+
+   input_file_x = []
+   input_file_y = []  
+   input_file_x_UTM = []
+   input_file_y_UTM = []  
+   input_file_lat = []
+   input_file_lon = []  
+   
+   input_file_zt = []
+   input_file_buildings_2d = []
+   input_file_bridges_2d = []
+   input_file_building_id = []
+   input_file_bridges_id = []   
+   input_file_building_type = []
+   input_file_vegetation_type = []
+   input_file_vegetation_height = []   
+   input_file_pavement_type = []
+   input_file_water_type = []
+   input_file_street_type = []  
+   input_file_street_crossings = []     
+   input_file_soil_type = []
+
+   
+# Load all user parameters from config file
+   for i in range(0,len(config.sections())):
+
+      read_tmp = config.sections()[i]
+      
+      if ( read_tmp == 'global' ):
+         global_acronym = config.get(read_tmp, 'acronym')
+         global_angle = float(config.get(read_tmp, 'rotation_angle'))
+         global_author = config.get(read_tmp, 'author')
+         global_campaign = config.get(read_tmp, 'campaign')
+         global_comment = config.get(read_tmp, 'comment')
+         global_contact = config.get(read_tmp, 'contact_person')
+         global_data_content = config.get(read_tmp, 'data_content')   
+         global_dependencies = config.get(read_tmp, 'dependencies')
+         global_institution = config.get(read_tmp, 'institution')
+         global_keywords = config.get(read_tmp, 'keywords')
+         global_location = config.get(read_tmp, 'location')
+         global_palm_version = float(config.get(read_tmp, 'palm_version'))
+         global_references = config.get(read_tmp, 'references')
+         global_site = config.get(read_tmp, 'site')
+         global_source = config.get(read_tmp, 'source')
+         global_version = float(config.get(read_tmp, 'version'))
+
+      if ( read_tmp == 'settings' ):
+         settings_filename_out = config.get(read_tmp, 'filename_out')
+         settings_lai_season = config.get(read_tmp, 'lai_season') 
+         settings_bridge_width = float(config.get(read_tmp, 'bridge_width'))      
+ 
+      if ( read_tmp.split("_")[0] == 'domain' ):
+         ndomains = ndomains + 1
+         domain_names.append(read_tmp.split("_")[1])
+         domain_px.append(float(config.get(read_tmp, 'pixel_size')))
+         domain_nx.append(int(config.get(read_tmp, 'nx')))
+         domain_ny.append(int(config.get(read_tmp, 'ny')))
+         domain_dz.append(float(config.get(read_tmp, 'dz')))       
+         domain_3d.append(config.getboolean(read_tmp, 'buildings_3d'))
+         domain_hv.append(config.getboolean(read_tmp, 'allow_high_vegetation'))
+         domain_cg.append(config.getboolean(read_tmp, 'generate_vegetation_patches'))       
+         domain_ip.append(config.getboolean(read_tmp, 'interpolate_terrain'))
+         domain_za.append(config.getboolean(read_tmp, 'use_palm_z_axis')) 
+         if domain_ip[ndomains-1] and not domain_za[ndomains-1]:
+            domain_za[ndomains-1] = True
+            print("+++ Overwrite user setting for use_palm_z_axis")
+         
+         domain_parent.append(config.get(read_tmp, 'domain_parent'))
+
+         domain_x0.append(int(floor(float(config.get(read_tmp, 'origin_x'))/float(config.get(read_tmp, 'pixel_size')))))
+         domain_y0.append(int(floor(float(config.get(read_tmp, 'origin_y'))/float(config.get(read_tmp, 'pixel_size')))))       
+         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'))))
+         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'))))
+
+      if ( read_tmp.split("_")[0] == 'input' ):
+         input_names.append(read_tmp.split("_")[1])
+         input_px.append(float(config.get(read_tmp, 'pixel_size')))
+         input_file_x.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_x'))
+         input_file_y.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_y')) 
+         input_file_lat.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_lat'))
+         input_file_lon.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_lon'))   
+         input_file_x_UTM.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_x_UTM'))
+         input_file_y_UTM.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_y_UTM'))  
+         input_file_zt.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_zt'))
+         input_file_buildings_2d.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_buildings_2d'))
+         input_file_bridges_2d.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_bridges_2d'))
+         input_file_building_id.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_building_id'))  
+         input_file_bridges_id.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_bridges_id')) 
+         input_file_building_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_building_type'))  
+         input_file_vegetation_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_vegetation_type')) 
+         input_file_vegetation_height.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_vegetation_height'))          
+         input_file_pavement_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_pavement_type'))    
+         input_file_water_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_water_type'))
+         input_file_street_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_street_type'))    
+         input_file_street_crossings.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_street_crossings'))  
+         #input_file_soil_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_soil_type'))    
+   return 0
+
+
+############################################################
+
+# Start of main program
+datatypes = {
+   "x": "f4",  
+   "y": "f4", 
+   "z": "f4", 
+   "lat": "f4",  
+   "lon": "f4",   
+   "E_UTM": "f4",  
+   "N_UTM": "f4", 
+   "zt": "f4",
+   "buildings_2d": "f4",
+   "buildings_3d": "b",
+   "bridges_2d": "f4",
+   "building_id": "i",
+   "bridges_id": "i",
+   "building_type": "b",
+   "nsurface_fraction": "i",
+   "vegetation_type": "b",
+   "vegetation_height": "f4",
+   "pavement_type": "b", 
+   "water_type": "b", 
+   "street_type": "b",  
+   "street_crossings": "b",     
+   "soil_type": "b",
+   "surface_fraction": "f4"
+   }
+
+fillvalues = {
+   "lat": float(-9999.0),
+   "lon": float(-9999.0),
+   "E_UTM": float(-9999.0),
+   "N_UTM": float(-9999.0),   
+   "zt": float(-9999.0),
+   "buildings_2d": float(-9999.0),
+   "buildings_3d": np.byte(-127),
+   "bridges_2d": float(-9999.0),
+   "building_id": int(-9999),
+   "bridges_id": int(-9999),
+   "building_type": np.byte(-127),
+   "nsurface_fraction": int(-9999),
+   "vegetation_type": np.byte(-127),
+   "vegetation_height": float(-9999.0),
+   "pavement_type": np.byte(-127),
+   "water_type": np.byte(-127),
+   "street_type": np.byte(-127), 
+   "street_crossings": np.byte(-127),   
+   "soil_type": np.byte(-127),
+   "surface_fraction": float(-9999.0)
+   }
+
+defaultvalues = {
+   "lat": float(-9999.0),
+   "lon": float(-9999.0),
+   "E_UTM": float(-9999.0),
+   "N_UTM": float(-9999.0),  
+   "zt": float(0.0),
+   "buildings_2d": float(0.0),
+   "buildings_3d": 0,   
+   "bridges_2d": float(0.0),
+   "building_id": int(0),
+   "bridges_id": int(0),
+   "building_type": 1,
+   "nsurface_fraction": int(-9999),
+   "vegetation_type": 3,
+   "vegetation_height": float(-9999.0),
+   "pavement_type": 1, 
+   "water_type": 1, 
+   "street_type": 1, 
+   "street_crossings": 0,    
+   "soil_type": 1,
+   "surface_fraction": float(0.0)
+   }
+  
+# Read configuration file and set parameters accordingly
+read_config_file()
+
+
+filename = []
+ii = []
+ii_parent = []
+# Define indices and filenames for all domains and create netCDF files
+for i in range(0,ndomains):
+
+#  Calculate indices and input files
+   ii.append(input_px.index(domain_px[i]))
+   filename.append(settings_filename_out + "_" + domain_names[i])
+   if domain_parent[i] is not None:
+      ii_parent.append(domain_names.index(domain_parent[i]))
+   else:
+      ii_parent.append(None)
+
+
+   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])
+   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])
+   lat = nc_read_from_file_2d(input_file_lat[ii[i]], "Band1", domain_x0[i], domain_x0[i], domain_y0[i], domain_y0[i])
+   lon = nc_read_from_file_2d(input_file_lon[ii[i]], "Band1", domain_x0[i], domain_x0[i], domain_y0[i], domain_y0[i])  
+   
+#  Create NetCDF output file and set global attributes
+   nc_create_file(filename[i])
+   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,global_version)
+
+   
+# Process terrain height
+for i in range(0,ndomains):
+
+#  Read and write terrain height (zt)
+   zt = nc_read_from_file_2d(input_file_zt[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])
+
+#  Final step: add zt array to the global array   
+   zt_all.append(zt)
+   del zt
+ 
+#  Calculate the global (all domains) minimum of the terrain height. This value will be substracted for all
+#  data sets
+zt_min = min(zt_all[0].flatten())  
+for i in range(0,ndomains):
+   zt_min = min(zt_min,min(zt_all[i].flatten()))
+   
+del zt_all[:]
+
+for i in range(0,ndomains):
+
+#  Read and write terrain height (zt)
+   zt = nc_read_from_file_2d(input_file_zt[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])
+   x = nc_read_from_file_1d(input_file_x[ii[i]], "x", domain_x0[i], domain_x1[i])
+   y = nc_read_from_file_1d(input_file_y[ii[i]], "y", domain_y0[i], domain_y1[i])
+
+   print( "Shift terrain height by -" + str(zt_min))
+   zt = zt - zt_min
+
+#  If necessary, interpolate parent domain terrain height on child domain grid and blend the two
+   if domain_ip[i]:
+      tmp_x0 = int( domain_x0[i] * domain_px[i] / domain_px[ii_parent[i]] ) - 1
+      tmp_y0 = int( domain_y0[i] * domain_px[i] / domain_px[ii_parent[i]] ) - 1
+      tmp_x1 = int( domain_x1[i] * domain_px[i] / domain_px[ii_parent[i]] ) + 1
+      tmp_y1 = int( domain_y1[i] * domain_px[i] / domain_px[ii_parent[i]] ) + 1
+                
+      tmp_x = nc_read_from_file_1d(input_file_x[ii_parent[i]], "x", tmp_x0, tmp_x1)
+      tmp_y = nc_read_from_file_1d(input_file_y[ii_parent[i]], "y", tmp_y0, tmp_y1)
+
+      zt_parent = nc_read_from_file_2d(input_file_zt[ii_parent[i]], 'Band1', tmp_x0, tmp_x1, tmp_y0, tmp_y1)
+
+      print( "Shift terrain height by -" + str(zt_min))
+      zt_parent = zt_parent - zt_min
+
+#     Interpolate array and bring to PALM grid of child domain
+      zt_ip = interpolate_2d(zt_parent,tmp_x,tmp_y,x,y)
+      zt_ip = bring_to_palm_grid(zt_ip,x,y,domain_dz[ii_parent[i]])
+
+#     Shift the child terrain height according to the parent mean terrain height       
+      zt = zt - np.mean(zt)  + np.mean(zt_ip)
+  
+  
+#     Blend over the parent and child terrain height within a radius of 50 px  
+      zt = blend_array_2d(zt,zt_ip,50)
+   
+#  Final step: add zt array to the global array   
+   zt_all.append(zt)
+   del zt
+ 
+
+# Read and shift x and y coordinates, shift terrain height according to its minimum value and write all data
+# to file  
+for i in range(0,ndomains):   
+#  Read horizontal grid variables from zt file and write them to output file
+   x = nc_read_from_file_1d(input_file_x[ii[i]], "x", domain_x0[i], domain_x1[i])
+   y = nc_read_from_file_1d(input_file_y[ii[i]], "y", domain_y0[i], domain_y1[i])  
+   x = x - min(x.flatten()) + domain_px[i]/2.0
+   y = y - min(y.flatten()) + domain_px[i]/2.0  
+   nc_write_dimension(filename[i], 'x', x, datatypes["x"])
+   nc_write_dimension(filename[i], 'y', y, datatypes["y"])  
+   nc_write_attribute(filename[i], 'x', 'long_name', 'x')
+   nc_write_attribute(filename[i], 'x', 'standard_name','projection_x_coordinate')   
+   nc_write_attribute(filename[i], 'x', 'units', 'm')
+   nc_write_attribute(filename[i], 'y', 'long_name', 'x')
+   nc_write_attribute(filename[i], 'y', 'standard_name', 'projection_y_coordinate')   
+   nc_write_attribute(filename[i], 'y', 'units', 'm')
+
+   lat = nc_read_from_file_2d(input_file_lat[ii[i]], "Band1", domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])
+   lon = nc_read_from_file_2d(input_file_lon[ii[i]], "Band1", domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])  
+
+   nc_write_to_file_2d(filename[i], 'lat', lat, datatypes["lat"],'y','x',fillvalues["lat"])
+   nc_write_to_file_2d(filename[i], 'lon', lon, datatypes["lon"],'y','x',fillvalues["lon"])  
+   
+   nc_write_attribute(filename[i], 'lat', 'long_name', 'latitude')
+   nc_write_attribute(filename[i], 'lat', 'standard_name','latitude')   
+   nc_write_attribute(filename[i], 'lat', 'units', 'degrees_north')
+ 
+   nc_write_attribute(filename[i], 'lon', 'long_name', 'longitude')
+   nc_write_attribute(filename[i], 'lon', 'standard_name','longitude')   
+   nc_write_attribute(filename[i], 'lon', 'units', 'degrees_east')  
+   
+   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])
+   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])  
+   
+   nc_write_to_file_2d(filename[i], 'E_UTM', x_UTM, datatypes["E_UTM"],'y','x',fillvalues["E_UTM"])
+   nc_write_to_file_2d(filename[i], 'N_UTM', y_UTM, datatypes["N_UTM"],'y','x',fillvalues["N_UTM"])  
+   
+   nc_write_attribute(filename[i], 'E_UTM', 'long_name', 'easting')
+   nc_write_attribute(filename[i], 'E_UTM', 'standard_name','projection_x_coorindate')   
+   nc_write_attribute(filename[i], 'E_UTM', 'units', 'm')
+ 
+   nc_write_attribute(filename[i], 'N_UTM', 'long_name', 'northing')
+   nc_write_attribute(filename[i], 'N_UTM', 'standard_name','projection_y_coorindate')   
+   nc_write_attribute(filename[i], 'N_UTM', 'units', 'm')  
+  
+   nc_write_crs(filename[i])
+   
+  
+  
+#  If necessary, bring terrain height to PALM's vertical grid. This is either forced by the user or implicitly
+#  by using interpolation for a child domain        
+   if domain_za[i]:
+      zt_all[i] = bring_to_palm_grid(zt_all[i],x,y,domain_dz[i])
+
+   nc_write_to_file_2d(filename[i], 'zt', zt_all[i], datatypes["zt"],'y','x',fillvalues["zt"])
+   nc_write_attribute(filename[i], 'zt', 'long_name', 'orography')
+   nc_write_attribute(filename[i], 'zt', 'units', 'm')
+   nc_write_attribute(filename[i], 'zt', 'res_orig', domain_px[i]) 
+   nc_write_attribute(filename[i], 'zt', 'coordinates', 'E_UTM N_UTM lon lat') 
+   nc_write_attribute(filename[i], 'zt', 'grid_mapping', 'E_UTM N_UTM lon lat')   
+
+del zt_all
+
+
+#  Process building height, id, and type
+for i in range(0,ndomains):  
+   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])
+   
+   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])  
+ 
+   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])  
+   building_type[building_type == 255] = fillvalues["building_type"]
+   building_type = np.where(building_type < 1,defaultvalues["building_type"],building_type)
+ 
+   check = check_arrays_2(buildings_2d,building_id,fillvalues["buildings_2d"],fillvalues["building_id"])
+   if not check:
+      buildings_2d = np.where(building_id != fillvalues["building_id"],buildings_2d,fillvalues["buildings_2d"])
+      building_id = np.where(buildings_2d == fillvalues["buildings_2d"],fillvalues["building_id"],building_id)
+      print("Data check #1 " + str(check_arrays_2(buildings_2d,building_id,fillvalues["buildings_2d"],fillvalues["building_id"])))
+
+   check = check_arrays_2(buildings_2d,building_type,fillvalues["buildings_2d"],fillvalues["building_type"])
+   if not check:
+      building_type = np.where(buildings_2d == fillvalues["buildings_2d"],fillvalues["building_type"],building_type)
+      building_type = np.where((building_type == fillvalues["building_type"]) & (buildings_2d != fillvalues["buildings_2d"]),defaultvalues["building_type"],building_type)
+      print("Data check #2 " + str(check_arrays_2(buildings_2d,building_type,fillvalues["buildings_2d"],fillvalues["building_type"])))
+ 
+   nc_write_to_file_2d(filename[i], 'buildings_2d', buildings_2d, datatypes["buildings_2d"],'y','x',fillvalues["buildings_2d"])
+   nc_write_attribute(filename[i], 'buildings_2d', 'long_name', 'buildings')
+   nc_write_attribute(filename[i], 'buildings_2d', 'units', 'm')
+   nc_write_attribute(filename[i], 'buildings_2d', 'res_orig', domain_px[i])  
+   nc_write_attribute(filename[i], 'buildings_2d', 'lod', 1) 
+   nc_write_attribute(filename[i], 'buildings_2d', 'coordinates', 'E_UTM N_UTM lon lat') 
+   nc_write_attribute(filename[i], 'buildings_2d', 'grid_mapping', 'E_UTM N_UTM lon lat')   
+   
+   nc_write_to_file_2d(filename[i], 'building_id', building_id, datatypes["building_id"],'y','x',fillvalues["building_id"])
+   nc_write_attribute(filename[i], 'building_id', 'long_name', 'building id')
+   nc_write_attribute(filename[i], 'building_id', 'units', '')
+   nc_write_attribute(filename[i], 'building_id', 'res_orig', domain_px[i]) 
+   nc_write_attribute(filename[i], 'building_id', 'coordinates', 'E_UTM N_UTM lon lat') 
+   nc_write_attribute(filename[i], 'building_id', 'grid_mapping', 'E_UTM N_UTM lon lat')   
+   
+   nc_write_to_file_2d(filename[i], 'building_type', building_type, datatypes["building_type"],'y','x',fillvalues["building_type"])
+   nc_write_attribute(filename[i], 'building_type', 'long_name', 'building type')
+   nc_write_attribute(filename[i], 'building_type', 'units', '')
+   nc_write_attribute(filename[i], 'building_type', 'res_orig', domain_px[i])  
+   nc_write_attribute(filename[i], 'building_type', 'coordinates', 'E_UTM N_UTM lon lat') 
+   nc_write_attribute(filename[i], 'building_type', 'grid_mapping', 'E_UTM N_UTM lon lat')   
+ 
+del buildings_2d
+del building_id
+del building_type
+
+# Create 3d buildings if necessary. In that course, read bridge objects and add them to building layer
+for i in range(0,ndomains): 
+   
+   if domain_3d[i]:
+      x = nc_read_from_file_2d_all(filename[i], 'x')
+      y = nc_read_from_file_2d_all(filename[i], 'y')  
+      buildings_2d = nc_read_from_file_2d_all(filename[i], 'buildings_2d') 
+      building_id = nc_read_from_file_2d_all(filename[i], 'building_id') 
+      
+      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])
+      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])
+      
+      bridges_2d = np.where(bridges_2d == 0.0,fillvalues["bridges_2d"],bridges_2d)
+      building_id = np.where(bridges_2d == fillvalues["bridges_2d"],building_id,bridges_id)
+      
+      
+      if np.any(buildings_2d != fillvalues["buildings_2d"]):
+         buildings_3d, z = make_3d_from_2d(buildings_2d,x,y,domain_dz[i])
+         if np.any(bridges_2d != fillvalues["bridges_2d"]):
+            buildings_3d = make_3d_from_bridges_2d(buildings_3d,bridges_2d,x,y,domain_dz[i],settings_bridge_width,fillvalues["bridges_2d"])
+         else:
+            print("Skipping creation of 3D bridges (no bridges in domain)")
+            
+            
+         nc_write_dimension(filename[i], 'z', z, datatypes["z"])  
+         nc_write_attribute(filename[i], 'z', 'long_name', 'z')  
+         nc_write_attribute(filename[i], 'z', 'units', 'm')
+ 
+         nc_write_to_file_3d(filename[i], 'buildings_3d', buildings_3d, datatypes["buildings_3d"],'z','y','x',fillvalues["buildings_3d"])    
+         nc_write_attribute(filename[i], 'buildings_3d', 'long_name', 'buildings 3d')
+         nc_write_attribute(filename[i], 'buildings_3d', 'units', '')
+         nc_write_attribute(filename[i], 'buildings_3d', 'res_orig', domain_px[i])  
+         nc_write_attribute(filename[i], 'buildings_3d', 'lod', 2)  
+         
+         del buildings_3d
+         
+      else:
+         print("Skipping creation of 3D buildings (no buildings in domain)")
+
+
+      del bridges_2d, bridges_id, building_id
+
+
+
+# Read vegetation type, water_type, pavement_type, soil_type and make fields consistent
+for i in range(0,ndomains):  
+
+   building_type = nc_read_from_file_2d_all(filename[i], 'building_type') 
+     
+   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])   
+   vegetation_type[vegetation_type == 255] = fillvalues["vegetation_type"]
+   vegetation_type = np.where((vegetation_type < 1) & (vegetation_type != fillvalues["vegetation_type"]),defaultvalues["vegetation_type"],vegetation_type)
+   
+   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])   
+   pavement_type[pavement_type == 255] = fillvalues["pavement_type"]
+   pavement_type = np.where((pavement_type < 1) & (pavement_type != fillvalues["pavement_type"]),defaultvalues["pavement_type"],pavement_type)
+
+   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])   
+   water_type[water_type == 255] = fillvalues["water_type"]
+   water_type = np.where((water_type < 1) & (water_type != fillvalues["water_type"]),defaultvalues["water_type"],water_type)
+  
+#  to do: replace by real soil input data
+   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])  
+   soil_type[soil_type == 255] = fillvalues["soil_type"]
+   soil_type = np.where((soil_type < 1) & (soil_type != fillvalues["soil_type"]),defaultvalues["soil_type"],soil_type)
+
+#  Make arrays consistent
+#  #1 Set vegetation type to missing for pixel where a pavement type is set
+   vegetation_type = np.where((vegetation_type != fillvalues["vegetation_type"]) & (pavement_type != fillvalues["pavement_type"]),fillvalues["vegetation_type"],vegetation_type)
+
+#  #2 Set vegetation type to missing for pixel where a building type is set
+   vegetation_type = np.where((vegetation_type != fillvalues["vegetation_type"]) & (building_type != fillvalues["building_type"]) ,fillvalues["vegetation_type"],vegetation_type)
+
+#  #3 Set vegetation type to missing for pixel where a building type is set
+   vegetation_type = np.where((vegetation_type != fillvalues["vegetation_type"]) & (water_type != fillvalues["water_type"]),fillvalues["vegetation_type"],vegetation_type)   
+
+#  #4 Remove pavement for pixels with buildings 
+   pavement_type = np.where((pavement_type != fillvalues["pavement_type"]) & (building_type != fillvalues["building_type"]),fillvalues["pavement_type"],pavement_type)  
+
+#  #5 Remove pavement for pixels with water
+   pavement_type = np.where((pavement_type != fillvalues["pavement_type"]) & (water_type != fillvalues["water_type"]),fillvalues["pavement_type"],pavement_type)   
+
+#  #6 Remove water for pixels with buildings
+   water_type = np.where((water_type != fillvalues["water_type"]) & (building_type != fillvalues["building_type"]),fillvalues["water_type"],water_type) 
+   
+
+#  #7 to be removed: set default soil type everywhere
+   soil_type = np.where((vegetation_type != fillvalues["vegetation_type"]) | (pavement_type != fillvalues["pavement_type"]),defaultvalues["soil_type"],fillvalues["soil_type"]) 
+  
+  
+#  Check for consistency and fill empty fields with default vegetation type
+   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"])
+   
+   if test:
+      vegetation_type = np.where(consistency_array == 0,defaultvalues["vegetation_type"],vegetation_type)
+      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"])
+   
+#  Create surface_fraction array
+   x = nc_read_from_file_2d_all(filename[i], 'x')
+   y = nc_read_from_file_2d_all(filename[i], 'y')  
+   nsurface_fraction = np.arange(0,3)
+   surface_fraction = np.ones((len(nsurface_fraction),len(y),len(x)))
+   
+   surface_fraction[0,:,:] = np.where(vegetation_type != fillvalues["vegetation_type"], 1.0, 0.0)
+   surface_fraction[1,:,:] = np.where(pavement_type != fillvalues["pavement_type"], 1.0, 0.0)
+   surface_fraction[2,:,:] = np.where(water_type != fillvalues["water_type"], 1.0, 0.0)
+   
+   nc_write_dimension(filename[i], 'nsurface_fraction', nsurface_fraction, datatypes["nsurface_fraction"]) 
+   nc_write_to_file_3d(filename[i], 'surface_fraction', surface_fraction, datatypes["surface_fraction"],'nsurface_fraction','y','x',fillvalues["surface_fraction"])     
+   nc_write_attribute(filename[i], 'surface_fraction', 'long_name', 'surface fraction')
+   nc_write_attribute(filename[i], 'surface_fraction', 'units', '')
+   nc_write_attribute(filename[i], 'surface_fraction', 'res_orig', domain_px[i])     
+   del surface_fraction   
+   
+ 
+#  Correct vegetation_type when a vegetation height is available and is indicative of low vegeetation
+   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])   
+         
+   vegetation_type = np.where((vegetation_height != fillvalues["vegetation_height"]) & (vegetation_height == 0.0) & ((vegetation_type == 4) | (vegetation_type == 7) | (vegetation_type == 17)), 3, vegetation_type)
+   vegetation_height = np.where((vegetation_height != fillvalues["vegetation_height"]) & (vegetation_height == 0.0) & ((vegetation_type == 4) | (vegetation_type == 7) | (vegetation_type == 17)), fillvalues["vegetation_height"],vegetation_height)
+
+
+   nc_write_to_file_2d(filename[i], 'vegetation_type', vegetation_type, datatypes["vegetation_type"],'y','x',fillvalues["vegetation_type"])     
+   nc_write_attribute(filename[i], 'vegetation_type', 'long_name', 'vegetation type')
+   nc_write_attribute(filename[i], 'vegetation_type', 'units', '')
+   nc_write_attribute(filename[i], 'vegetation_type', 'res_orig', domain_px[i])  
+   nc_write_attribute(filename[i], 'vegetation_type', 'coordinates', 'E_UTM N_UTM lon lat') 
+   nc_write_attribute(filename[i], 'vegetation_type', 'grid_mapping', 'E_UTM N_UTM lon lat')   
+   del vegetation_type
+
+   nc_write_to_file_2d(filename[i], 'pavement_type', pavement_type, datatypes["pavement_type"],'y','x',fillvalues["pavement_type"])    
+   nc_write_attribute(filename[i], 'pavement_type', 'long_name', 'pavement type')
+   nc_write_attribute(filename[i], 'pavement_type', 'units', '')
+   nc_write_attribute(filename[i], 'pavement_type', 'res_orig', domain_px[i])  
+   nc_write_attribute(filename[i], 'pavement_type', 'coordinates', 'E_UTM N_UTM lon lat') 
+   nc_write_attribute(filename[i], 'pavement_type', 'grid_mapping', 'E_UTM N_UTM lon lat')   
+   del pavement_type
+
+   nc_write_to_file_2d(filename[i], 'water_type', water_type, datatypes["water_type"],'y','x',fillvalues["water_type"]) 
+   nc_write_attribute(filename[i], 'water_type', 'long_name', 'water type')
+   nc_write_attribute(filename[i], 'water_type', 'units', '')
+   nc_write_attribute(filename[i], 'water_type', 'res_orig', domain_px[i]) 
+   nc_write_attribute(filename[i], 'water_type', 'coordinates', 'E_UTM N_UTM lon lat') 
+   nc_write_attribute(filename[i], 'water_type', 'grid_mapping', 'E_UTM N_UTM lon lat')   
+   del water_type
+
+   nc_write_to_file_2d(filename[i], 'soil_type', soil_type, datatypes["soil_type"],'y','x',fillvalues["soil_type"])  
+   nc_write_attribute(filename[i], 'soil_type', 'long_name', 'soil type')
+   nc_write_attribute(filename[i], 'soil_type', 'units', '')
+   nc_write_attribute(filename[i], 'soil_type', 'res_orig', domain_px[i])   
+   nc_write_attribute(filename[i], 'soil_type', 'coordinates', 'E_UTM N_UTM lon lat') 
+   nc_write_attribute(filename[i], 'soil_type', 'grid_mapping', 'E_UTM N_UTM lon lat')   
+   del soil_type
+
+
+
+
+
+#  pixels with bridges get building_type = 7 = bridge. This does not change the _type setting for the under-bridge
+#  area
+   if domain_3d[i]:
+      if np.any(building_type != fillvalues["building_type"]):
+ 
+         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])         
+         bridges_2d = np.where(bridges_2d == 0.0,fillvalues["bridges_2d"],bridges_2d)
+         building_type = np.where(bridges_2d != fillvalues["bridges_2d"],7,building_type)
+         nc_overwrite_to_file_2d(filename[i], 'building_type', building_type)  
+ 
+         del building_type  
+         del bridges_2d
+
+# Read/Write street type and street crossings
+for i in range(0,ndomains):  
+
+   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])   
+   street_type[street_type == 255] = fillvalues["street_type"]
+   street_type = np.where((street_type < 1) & (street_type != fillvalues["street_type"]),defaultvalues["street_type"],street_type)
+   
+   nc_write_to_file_2d(filename[i], 'street_type', street_type, datatypes["street_type"],'y','x',fillvalues["street_type"]) 
+   nc_write_attribute(filename[i], 'street_type', 'long_name', 'street type')
+   nc_write_attribute(filename[i], 'street_type', 'units', '')
+   nc_write_attribute(filename[i], 'street_type', 'res_orig', domain_px[i])  
+   nc_write_attribute(filename[i], 'street_type', 'coordinates', 'E_UTM N_UTM lon lat') 
+   nc_write_attribute(filename[i], 'street_type', 'grid_mapping', 'E_UTM N_UTM lon lat')   
+   del street_type
+   
+   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])   
+   street_crossings[street_crossings == 255] = fillvalues["street_crossings"]
+   street_crossings = np.where((street_crossings < 1) & (street_crossings != fillvalues["street_crossings"]),defaultvalues["street_crossings"],street_crossings)
+   
+   nc_write_to_file_2d(filename[i], 'street_crossings', street_crossings, datatypes["street_crossings"],'y','x',fillvalues["street_crossings"]) 
+   nc_write_attribute(filename[i], 'street_crossings', 'long_name', 'street crossings')
+   nc_write_attribute(filename[i], 'street_crossings', 'units', '')
+   nc_write_attribute(filename[i], 'street_crossings', 'res_orig', domain_px[i])  
+   nc_write_attribute(filename[i], 'street_crossings', 'coordinates', 'E_UTM N_UTM lon lat') 
+   nc_write_attribute(filename[i], 'street_crossings', 'grid_mapping', 'E_UTM N_UTM lon lat')   
+   del street_crossings
Index: /palm/trunk/SCRIPTS/palm_csd_files/palm_csd_netcdf_interface.py
===================================================================
--- /palm/trunk/SCRIPTS/palm_csd_files/palm_csd_netcdf_interface.py	(revision 3567)
+++ /palm/trunk/SCRIPTS/palm_csd_files/palm_csd_netcdf_interface.py	(revision 3567)
@@ -0,0 +1,272 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+#--------------------------------------------------------------------------------#
+# This file is part of the PALM model system.
+#
+# PALM is free software: you can redistribute it and/or modify it under the terms
+# of the GNU General Public License as published by the Free Software Foundation,
+# either version 3 of the License, or (at your option) any later version.
+#
+# PALM is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+# A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along with
+# PALM. If not, see <http://www.gnu.org/licenses/>.
+#
+# Copyright 1997-2018  Leibniz Universitaet Hannover
+#--------------------------------------------------------------------------------#
+#
+# Current revisions:
+# -----------------
+# 
+# 
+# Former revisions:
+# -----------------
+# $Id$
+# Initial revisions
+#
+#
+# 
+#
+#
+# Description:
+# ------------
+# NetCDF interface routines for palm_csd
+#
+# @Author Bjoern Maronga (maronga@muk.uni-hannover.de)
+#
+#------------------------------------------------------------------------------#
+
+from netCDF4 import Dataset
+
+def nc_read_from_file_2d(filename, varname, x0, x1, y0, y1):
+
+
+   import numpy as np
+   import sys
+   
+   try:
+      f = open(filename)
+      f.close()
+#      print("Load: " + filename + ".")
+   except FileNotFoundError:
+      print("Error: " + filename + ". No such file. Aborting...")
+      sys.exit(1)
+   
+   nc_file = Dataset(filename, "r+", format="NETCDF4")
+   tmp_array = np.array(nc_file.variables[varname][y0:y1+1,x0:x1+1] , dtype=type(nc_file.variables[varname]))
+
+   return tmp_array
+
+def nc_read_from_file_2d_all(filename, varname):
+
+
+   import numpy as np
+   import sys
+   
+   try:
+      f = open(filename)
+      f.close()
+#      print("Load: " + filename + ".")
+   except FileNotFoundError:
+      print("Error: " + filename + ". No such file. Aborting...")
+      sys.exit(1)
+   
+   nc_file = Dataset(filename, "r+", format="NETCDF4")
+   tmp_array = np.array(nc_file.variables[varname][:][:], dtype=type(nc_file.variables[varname]))
+
+   return tmp_array
+
+def nc_read_from_file_1d(filename, varname, x0, x1):
+
+   import numpy as np
+   import sys
+   
+   try:
+      f = open(filename)
+      f.close()
+      #print("Load: " + filename + ".")
+   except FileNotFoundError:
+      print("Error: " + filename + ". No such file. Aborting...")
+      sys.exit(1)
+   
+   nc_file = Dataset(filename, "r+", format="NETCDF4")
+   tmp_array = np.array(nc_file.variables[varname][x0:x1+1] , dtype=type(nc_file.variables[varname]))
+
+   return tmp_array
+
+
+def nc_create_file(filename):
+
+   import datetime
+
+   try:
+      f =  Dataset(filename, "w", format="NETCDF4")
+      f.close()
+      print("Created: " + filename + ".")
+   except FileNotFoundError:
+      print("Error. Could not create file: " + filename + ". Aborting...")
+      sys.exit(1)
+
+   return 0
+
+
+   
+def nc_write_global_attributes(filename,origin_x,origin_y,origin_lat,origin_lon,origin_time, 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,global_version):
+
+   import datetime
+   import os
+
+   print("Writing global attributes to file...")
+
+   f =  Dataset(filename, "a", format="NETCDF4")
+
+   f.setncattr('Conventions',"CF-1.7") 
+   
+   f.origin_x = origin_x
+   f.origin_y = origin_y
+   f.origin_time = origin_time
+   f.origin_lat = origin_lat
+   f.origin_lon = origin_lon
+   
+   f.acronym = global_acronym
+   f.rotation_angle = global_angle
+   f.author = global_author
+   f.campaign = global_campaign
+   f.comment = global_comment
+   f.contact = global_contact
+   f.data_content   = global_data_content
+   f.dependencies = global_dependencies
+   f.institution = global_institution
+   f.keywords = global_keywords
+   f.location = global_location
+   f.palm_version = global_palm_version
+   f.references = global_references
+   f.site = global_site
+   f.source = global_source
+   f.version = global_version
+   
+
+   f.close()
+
+   return 0
+
+
+def nc_write_dimension(filename,varname,array,datatype):
+
+   try:
+      f =  Dataset(filename, "a", format="NETCDF4")
+      #print("Opened: " + filename + ".")
+   except FileNotFoundError:
+      print("Error. Could not open file: " + filename + ". Aborting...")
+      sys.exit(1)
+
+   print("Writing dimension " + varname + " to file...")
+   
+   f.createDimension(varname,len(array))
+   temp = f.createVariable(varname,datatype,varname)
+   temp[:] = array
+   
+   f.close()
+   
+   return 0
+
+def nc_write_to_file_2d(filename,varname,array,datatype,dimname1,dimname2,fillvalue):
+
+   try:
+      f =  Dataset(filename, "a", format="NETCDF4")
+      #print("Opened: " + filename + ".")
+   except FileNotFoundError:
+      print("Error. Could not open file: " + filename + ". Aborting...")
+      sys.exit(1)
+
+   print("Writing array " + varname + " to file...")
+   
+   temp = f.createVariable(varname,datatype,(dimname1,dimname2),fill_value=fillvalue)
+   temp[:] = array
+   
+   f.close()
+
+   return 0
+
+def nc_overwrite_to_file_2d(filename,varname,array):
+
+   try:
+      f =  Dataset(filename, "a", format="NETCDF4")
+      #print("Opened: " + filename + ".")
+   except FileNotFoundError:
+      print("Error. Could not open file: " + filename + ". Aborting...")
+      sys.exit(1)
+
+   print("Writing array " + varname + " to file...")
+   
+   temp = f.variables[varname]
+   temp[:,:] = array
+   
+   f.close()
+
+   return 0
+
+def nc_write_to_file_3d(filename,varname,array,datatype,dimname1,dimname2,dimname3,fillvalue):
+
+   try:
+      f =  Dataset(filename, "a", format="NETCDF4")
+      #print("Opened: " + filename + ".")
+   except FileNotFoundError:
+      print("Error. Could not open file: " + filename + ". Aborting...")
+      sys.exit(1)
+
+   print("Writing array " + varname + " to file...")
+   
+   temp = f.createVariable(varname,datatype,(dimname1,dimname2,dimname3),fill_value=fillvalue)
+   temp[:,:,:] = array
+   
+   f.close()
+
+   return 0
+
+
+def nc_write_attribute(filename,variable,attribute,value):
+
+   f = Dataset(filename, "a", format="NETCDF4")
+
+   var = f.variables[variable]
+   var.setncattr(attribute,value) 
+   
+   f.close()
+
+   return 0
+
+def nc_write_crs(filename):
+
+   try:
+      f =  Dataset(filename, "a", format="NETCDF4")
+      #print("Opened: " + filename + ".")
+   except FileNotFoundError:
+      print("Error. Could not open file: " + filename + ". Aborting...")
+      sys.exit(1)
+
+   print("Writing crs to file...")
+   
+   temp = f.createVariable("crs","i")
+   
+   temp.long_name = "coordinate reference system"
+   temp.grid_mapping_name = "transverse_mercator"
+   temp.semi_major_axis = 6378137.0
+   temp.inverse_flattening = 298.257222101
+   temp.longitude_of_prime_meridian = 0.0
+   temp.longitude_of_central_meridian = 9.0
+   temp.scale_factor_at_central_meridian = 0.9996
+   temp.latitude_of_projection_origin = 0.0
+   temp.false_easting = 50000.0
+   temp.false_northing = 0.0
+   temp.units = "m"
+   temp.epsg_code = "EPSG:25833"
+   
+   f.close()
+
+   return 0
Index: /palm/trunk/SCRIPTS/palm_csd_files/palm_csd_tools.py
===================================================================
--- /palm/trunk/SCRIPTS/palm_csd_files/palm_csd_tools.py	(revision 3567)
+++ /palm/trunk/SCRIPTS/palm_csd_files/palm_csd_tools.py	(revision 3567)
@@ -0,0 +1,150 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+#--------------------------------------------------------------------------------#
+# This file is part of the PALM model system.
+#
+# PALM is free software: you can redistribute it and/or modify it under the terms
+# of the GNU General Public License as published by the Free Software Foundation,
+# either version 3 of the License, or (at your option) any later version.
+#
+# PALM is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+# A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along with
+# PALM. If not, see <http://www.gnu.org/licenses/>.
+#
+# Copyright 1997-2018  Leibniz Universitaet Hannover
+#--------------------------------------------------------------------------------#
+#
+# Current revisions:
+# -----------------
+# 
+# 
+# Former revisions:
+# -----------------
+# $Id$
+# Initial revisions
+#
+#
+# 
+#
+#
+# Description:
+# ------------
+# Support routines for palm_csd
+#
+# @Author Bjoern Maronga (maronga@muk.uni-hannover.de)
+#
+#------------------------------------------------------------------------------#
+import numpy as np
+from scipy.interpolate import interp2d
+
+def blend_array_2d(array1,array2,radius):
+#  Blend over the parent and child terrain height within a radius of 50 px
+   
+   gradient_matrix = np.copy(array1)
+   gradient_matrix[:,:] = 1.0
+   gradient_px     = 50
+
+   for i in range(0,gradient_px):
+      gradient_matrix[:,i] = float(i)/float(gradient_px)
+
+
+   for i in range(len(gradient_matrix[0,:])-gradient_px,len(gradient_matrix[0,:])):
+      gradient_matrix[:,i] = float(len(gradient_matrix[0,:])-i)/float(gradient_px)  
+     
+
+   for j in range(0,gradient_px):
+      for i in range(0,len(gradient_matrix[0,:])):
+         if  gradient_matrix[j,i] == 1.0:
+            gradient_matrix[j,i] = float(j)/float(gradient_px) 
+         else:
+            gradient_matrix[j,i] = (gradient_matrix[j,i] + float(j)/float(gradient_px))/2.0
+
+
+   for j in range(len(gradient_matrix[:,0])-gradient_px,len(gradient_matrix[:,0])):
+      for i in range(0,len(gradient_matrix[0,:])):        
+         if  gradient_matrix[j,i] == 1.0:
+            gradient_matrix[j,i] = (len(gradient_matrix[:,0])-j)/float(gradient_px)
+         else:
+            gradient_matrix[j,i] = (gradient_matrix[j,i] + (len(gradient_matrix[:,0])-j)/float(gradient_px))/2.0
+
+   array_blended = array1 * gradient_matrix + (1.0 - gradient_matrix ) * array2 
+       
+   return array_blended
+
+
+def interpolate_2d(array,x1,y1,x2,y2):
+    
+      tmp_int2d = interp2d(y1,x1,array,kind='linear')
+      array_ip = tmp_int2d(y2.astype(float), x2.astype(float)) 
+
+      return array_ip 
+   
+   
+def bring_to_palm_grid(array,x,y,dz):
+    
+#     Bring the parent terrain height to the child grid  
+      k_tmp = np.arange(0,max(array.flatten())+dz*2,dz)
+      k_tmp[1:] = k_tmp[1:] - dz * 0.5
+      for l in range(0,len(x)):
+         for m in range(0,len(y)):
+            for k in range(0,len(k_tmp+1)):
+               if k_tmp[k] > array[m,l]:
+                  array[m,l] = k_tmp[k]-0.5*dz
+                  break
+
+      return array  
+
+
+def make_3d_from_2d(array_2d,x,y,dz):
+    
+      k_tmp = np.arange(0,max(array_2d.flatten())+dz*2,dz)
+ 
+      k_tmp[1:] = k_tmp[1:] - dz * 0.5
+      array_3d = np.ones((len(k_tmp),len(y),len(x)))
+      
+      for l in range(0,len(x)-1):
+         for m in range(0,len(y)-1):
+            for k in range(0,len(k_tmp)-1):
+               if k_tmp[k] > array_2d[m,l]:
+                  array_3d[k,m,l] = 0
+
+      return array_3d.astype(np.byte), k_tmp
+
+
+def make_3d_from_bridges_2d(array_3d,array_2d,x,y,dz,width,fill):
+      
+      for l in range(0,len(x)-1):
+         for m in range(0,len(y)-1):
+            if array_2d[m,l] != fill:
+               k_min = max( int(array_2d[m,l] - width)/dz, 0 )
+               k_max = int(round(array_2d[m,l]/dz))
+               array_3d[k_min:k_max+1,m,l] = 1
+
+
+      return array_3d.astype(np.byte)
+
+
+def check_arrays_2(array1,array2,fill1,fill2):
+
+   missing1 = np.where(array1 == fill1,1,0)
+   missing2 = np.where(array2 == fill2,1,0)
+   result = np.array_equal(missing1,missing2)
+    
+   return result
+
+
+def check_consistency_4(array1,array2,array3,array4,fill1,fill2,fill3,fill4):
+
+   tmp_array = np.where(array1 != fill1,1,0) + np.where(array2 != fill2,1,0) + np.where(array3 != fill3,1,0) + np.where(array4 != fill4,1,0)
+   
+   test = np.any(tmp_array.flatten() != 1)
+   if test:
+      print("*_type arrays are not consistent!")
+      print("max: " + str(max(tmp_array.flatten())) + ", min: " + str(min(tmp_array.flatten())))
+   else:
+      print("*_type arrays are consistent!")     
+   return tmp_array, test
+
Index: /palm/trunk/SCRIPTS/palm_wd
===================================================================
--- /palm/trunk/SCRIPTS/palm_wd	(revision 3566)
+++ /palm/trunk/SCRIPTS/palm_wd	(revision 3567)
@@ -25,4 +25,8 @@
 # -----------------
 # $Id$
+# Adjustments and bugfixes 
+#
+# 
+# 2825 2018-02-20 21:48:27Z maronga
 # Modified header
 # 
@@ -91,27 +95,24 @@
 
 # Determine PALM directories
-try: 
-   devnull = open(os.devnull, 'w')
-   out = sub.check_output("echo $PALM_BIN", shell=True, stderr=sub.STDOUT)
-   palm_bin = out.rstrip()
-   palm_dir = out.split("palm")[0] + "palm/" + out.split("palm")[1].split("/")[1]
-   out = None
-except:   
-   print "Error. $PALM_BIN is not set."
-   raise SystemExit
-
-
-# Read configuration file
-# First check if the configuration file exists
-if ( os.path.exists(palm_dir + '/.wd.config') == False ):
-    print "Error. No configuration file .wd.config found in " + palm_dir
-    raise SystemExit     
-
-config = ConfigParser.RawConfigParser()
-config.read(palm_dir + '/.wd.config')
+#  Check if configuration files exists and quit otherwise
+input_config = ".wd.config"
+for i in range(1,len(sys.argv)): 
+   input_config = str(sys.argv[i])
+        
+config = configparser.RawConfigParser(allow_no_value=True)
+
+if ( os.path.isfile(input_config) == False ):
+   print ("Error. No configuration file " + input_config + " found.")
+   raise SystemExit   
+else:
+   print(os.path.isfile(input_config))
+
+config.read(input_config)  
+
 
 description = []
 hostname = []
 username = []
+pycall = []
 
 for i in range(0,len(config.sections())):
@@ -125,4 +126,5 @@
           hostname.append(config.get(description_tmp, 'hostname'))
           username.append(config.get(description_tmp, 'username')) 
+          pycall.append(config.get(description_tmp, 'pycall'))
     else:
        update_frequency = int(config.get(description_tmp, 'update_frequency'))*60000
@@ -401,5 +403,5 @@
 	    # Perform ssh command
             host    = username[h] + "@" + hostname[h]
-            ssh = sub.Popen(["ssh", "%s" % host, "/sw/tools/python/2.7.6/generic/bin/python palm_wdd queue " + username[h]],
+            ssh = sub.Popen(["ssh", "%s" % host, pycall[h] + " palm_wdd queue " + username[h]],
                            shell=False,
                            stdout=sub.PIPE,
@@ -692,5 +694,5 @@
                   host = username[h] + "@" + hostname[h]
 
-           ssh = sub.Popen(["ssh", "%s" % host, "/sw/tools/python/2.7.6/generic/bin/python palm_wdd cancel " + jobname],
+           ssh = sub.Popen(["ssh", "%s" % host, pycall[h] + " palm_wdd cancel " + jobname],
                               shell=False,
                               stdout=sub.PIPE,
@@ -719,5 +721,5 @@
                 host = username[h] + "@" + hostname[h]
 
-        ssh = sub.Popen(["ssh", "%s" % host, "/sw/tools/python/2.7.6/generic/bin/python palm_wdd check " + jobname],
+        ssh = sub.Popen(["ssh", "%s" % host, pycall[h] + " palm_wdd check " + jobname],
                            shell=False,
                            stdout=sub.PIPE,
@@ -746,5 +748,5 @@
                 user = username[h]
 
-        ssh = sub.Popen(["ssh", "%s" % host, "/sw/tools/python/2.7.6/generic/bin/python palm_wdd stop " + jobrealname],
+        ssh = sub.Popen(["ssh", "%s" % host, pycall[h] + " palm_wdd stop " + jobrealname],
                            shell=False,
                            stdout=sub.PIPE,
@@ -778,5 +780,5 @@
                 user = username[h]
 
-        ssh = sub.Popen(["ssh", "%s" % host, "/sw/tools/python/2.7.6/generic/bin/python palm_wdd restart " + jobrealname],
+        ssh = sub.Popen(["ssh", "%s" % host, pycall[h] + " palm_wdd restart " + jobrealname],
                            shell=False,
                            stdout=sub.PIPE,
@@ -810,5 +812,5 @@
                 user = username[h]
 
-        ssh = sub.Popen(["ssh", "%s" % host, "/sw/tools/python/2.7.6/generic/bin/python palm_wdd parin " + jobrealname],
+        ssh = sub.Popen(["ssh", "%s" % host, pycall[h] + " palm_wdd parin " + jobrealname],
                            shell=False,
                            stdout=sub.PIPE,
@@ -843,5 +845,5 @@
                 user = username[h]
 
-        ssh = sub.Popen(["ssh", "%s" % host, "/sw/tools/python/2.7.6/generic/bin/python palm_wdd rc " + jobrealname],
+        ssh = sub.Popen(["ssh", "%s" % host, pycall[h] + " palm_wdd rc " + jobrealname],
                            shell=False,
                            stdout=sub.PIPE,
