Home

Context Navigation

source: palm/trunk/SCRIPTS/palm_csd @ 3689

Last change on this file since 3689 was 3688, checked in by maronga, 6 years ago
bugfixes in palm_csd
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 54.9 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |