#!/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 . # # Copyright 1997-2018 Leibniz Universitaet Hannover #--------------------------------------------------------------------------------# # # Current revisions: # ----------------- # # # Former revisions: # ----------------- # $Id: palm_csd 3688 2019-01-22 10:44:20Z maronga $ # Some unspecified bugfixes # # 3668 2019-01-14 12:49:24Z maronga # Various improvements and bugfixes # # 3629 2018-12-13 12:18:54Z maronga # Added canopy generator calls. Some improvements # # 3567 2018-11-27 13:59:21Z maronga # Initial revisions # # Description: # ------------ # Processing tool for creating PIDS conform static drivers from rastered NetCDF # input # # @Author Bjoern Maronga (maronga@muk.uni-hannover.de) # # @todo Make input files optional # @todo Allow for ASCII input of terrain height and building height # @todo Modularize reading config file # @todo Convert to object-oriented treatment (buidings, trees) # @todo Automatically shift child domains so that their origin lies intersects # a edge note of the parent grid #------------------------------------------------------------------------------# from palm_csd_files.palm_csd_netcdf_interface import * from palm_csd_files.palm_csd_tools import * from palm_csd_files.palm_csd_canopy_generator 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 config.read(input_config) # Definition of settings global settings_bridge_width global settings_lai_roof_intensive global settings_lai_roof_extensive global settings_season global settings_lai_low_default global settings_lai_high_default global settings_patch_height_default global settings_lai_alpha global settings_lai_beta 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 path_out global filename_out global version_out # 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_high_vegetation global domain_ip global domain_za global domain_parent global domain_green_roofs global domain_street_trees global domain_canopy_patches # 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_lai 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 input_file_vegetation_on_roofs global input_file_tree_crown_diameter global input_file_tree_height global input_file_tree_trunk_diameter global input_file_tree_type global input_file_patch_height global zt_all global zt_min settings_bridge_width = 3.0 settings_lai_roof_intensive = 0.5 settings_lai_roof_extensive = 1.0 settings_season = "summer" settings_lai_high_default = 6.0 settings_lai_low_default = 1.0 settings_patch_height_default = 10.0 settings_lai_alpha = 5.0 settings_lai_beta = 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 = "" path_out = "" version_out = 1 filename_out = "default" domain_names = [] domain_px = [] domain_x0 = [] domain_y0 = [] domain_x1 = [] domain_y1 = [] domain_nx = [] domain_ny = [] domain_dz = [] domain_3d = [] domain_high_vegetation = [] domain_ip = [] domain_za = [] domain_parent = [] domain_green_roofs = [] domain_street_trees = [] domain_canopy_patches = [] 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_lai = [] 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 = [] input_file_vegetation_on_roofs = [] input_file_tree_crown_diameter = [] input_file_tree_height = [] input_file_tree_trunk_diameter = [] input_file_tree_type = [] input_file_patch_height = [] # 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') if ( read_tmp == 'settings' ): settings_lai_roof_intensive = config.get(read_tmp, 'lai_roof_intensive') settings_lai_roof_extensive = config.get(read_tmp, 'lai_roof_extensive') settings_bridge_width = float(config.get(read_tmp, 'bridge_width')) settings_season = config.get(read_tmp, 'season') settings_lai_high_default = float(config.get(read_tmp, 'lai_high_vegetation_default')) settings_lai_low_default = float(config.get(read_tmp, 'lai_low_vegetation_default')) settings_patch_height_default = float(config.get(read_tmp, 'patch_height_default')) settings_lai_alpha = float(config.get(read_tmp, 'lai_alpha')) settings_lai_beta = float(config.get(read_tmp, 'lai_beta')) if ( read_tmp == 'output' ): path_out = config.get(read_tmp, 'path') filename_out = config.get(read_tmp, 'file_out') version_out = float(config.get(read_tmp, 'version')) 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_high_vegetation.append(config.getboolean(read_tmp, 'allow_high_vegetation')) domain_canopy_patches.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')))) domain_green_roofs.append(config.getboolean(read_tmp, 'vegetation_on_roofs')) domain_street_trees.append(config.getboolean(read_tmp, 'street_trees')) 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_lai.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_lai')) 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_patch_height.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_patch_height')) tmp = config.get(read_tmp, 'file_tree_crown_diameter') if tmp is not None: input_file_tree_crown_diameter.append(config.get(read_tmp, 'path') + "/" + tmp) else: input_file_tree_crown_diameter.append(None) input_file_tree_height.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_tree_height')) input_file_tree_trunk_diameter.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_tree_trunk_diameter')) input_file_tree_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_tree_type')) input_file_vegetation_on_roofs.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_vegetation_on_roofs')) #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", "building_pars": "f4", "vegetation_pars": "f4", "tree_data": "f4", "tree_type": "b", "nbuilding_pars": "i", "nvegetation_pars": "i", "zlad": "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), "building_pars": float(-9999.0), "vegetation_pars": float(-9999.0), "tree_data": float(-9999.0), "tree_type": np.byte(-127) } 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), "buildings_pars": float(-9999.0), "tree_data": float(-9999.0), "tree_type": 0, "vegetation_pars": float(-9999.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(path_out + "/" + 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]) - 0.5 * domain_px[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]) - 0.5 * domain_px[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,version_out) del x_UTM, y_UTM, lat, lon # 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[:] print( "Shift terrain height by -" + str(zt_min)) 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]) zt = zt - zt_min nc_write_global_attribute(filename[i], 'origin_z', float(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 >= 254] = 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, buildings_2d # 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) # 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 == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18)), 3, vegetation_type) 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) # 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"]) # #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_3(vegetation_type,pavement_type,soil_type,fillvalues["vegetation_type"],fillvalues["pavement_type"],fillvalues["soil_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 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 del x del y # pixels with bridges get building_type = 7 = bridge. This does not change the _type setting for the under-bridge # area NOTE: when bridges are present the consistency check will fail at the moment 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 # Read/write vegetation on roofs for i in range(0,ndomains): if domain_green_roofs[i]: 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]) buildings_2d = nc_read_from_file_2d_all(filename[i], 'buildings_2d') x = nc_read_from_file_2d_all(filename[i], 'x') y = nc_read_from_file_2d_all(filename[i], 'y') nbuilding_pars = np.arange(0,46) building_pars = np.ones((len(nbuilding_pars),len(y),len(x))) building_pars[:,:,:] = fillvalues["building_pars"] # assign green fraction on roofs building_pars[3,:,:] = np.where( (buildings_2d != fillvalues["buildings_2d"] ) & ( green_roofs != fillvalues["building_pars"] ),1.0,fillvalues["building_pars"]) # assign leaf area index for vegetation on roofs building_pars[4,:,:] = np.where( ( buildings_2d != fillvalues["buildings_2d"] ) & ( green_roofs == 1.0 ),settings_lai_roof_intensive,fillvalues["building_pars"]) building_pars[4,:,:] = np.where( ( buildings_2d != fillvalues["buildings_2d"] ) & ( green_roofs == 2.0 ),settings_lai_roof_extensive,building_pars[4,:,:]) nc_write_dimension(filename[i], 'nbuilding_pars', nbuilding_pars, datatypes["nbuilding_pars"]) nc_write_to_file_3d(filename[i], 'building_pars', building_pars, datatypes["building_pars"],'nbuilding_pars','y','x',fillvalues["building_pars"]) nc_write_attribute(filename[i], 'building_pars', 'long_name', 'building_pars') nc_write_attribute(filename[i], 'building_pars', 'units', '') nc_write_attribute(filename[i], 'building_pars', 'res_orig', domain_px[i]) nc_write_attribute(filename[i], 'building_pars', 'coordinates', 'E_UTM N_UTM lon lat') nc_write_attribute(filename[i], 'building_pars', 'grid_mapping', 'E_UTM N_UTM lon lat') del building_pars, buildings_2d, x, y # Read tree data and create LAD and BAD arrays using the canopy generator for i in range(0,ndomains): lai = nc_read_from_file_2d(input_file_lai[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) vegetation_type = nc_read_from_file_2d_all(filename[i], 'vegetation_type') lai = np.where(vegetation_type == fillvalues["vegetation_type"],fillvalues["vegetation_pars"],lai) x = nc_read_from_file_2d_all(filename[i], 'x') y = nc_read_from_file_2d_all(filename[i], 'y') nvegetation_pars = np.arange(0,12) vegetation_pars = np.ones((len(nvegetation_pars),len(y),len(x))) vegetation_pars[:,:,:] = fillvalues["vegetation_pars"] vegetation_pars[1,:,:] = lai # Write out first version of LAI. Will later be overwritten. nc_write_dimension(filename[i], 'nvegetation_pars', nvegetation_pars, datatypes["nvegetation_pars"]) nc_write_to_file_3d(filename[i], 'vegetation_pars', vegetation_pars, datatypes["vegetation_pars"],'nvegetation_pars','y','x',fillvalues["vegetation_pars"]) nc_write_attribute(filename[i], 'vegetation_pars', 'long_name', 'vegetation_pars') nc_write_attribute(filename[i], 'vegetation_pars', 'units', '') nc_write_attribute(filename[i], 'vegetation_pars', 'res_orig', domain_px[i]) nc_write_attribute(filename[i], 'vegetation_pars', 'coordinates', 'E_UTM N_UTM lon lat') nc_write_attribute(filename[i], 'vegetation_pars', 'grid_mapping', 'E_UTM N_UTM lon lat') del lai, vegetation_pars, vegetation_type # Read tree data and create LAD and BAD arrays using the canopy generator for i in range(0,ndomains): if domain_street_trees[i]: vegetation_pars = nc_read_from_file_2d_all(filename[i], 'vegetation_pars') lai = np.copy(vegetation_pars[1,:,:]) x = nc_read_from_file_2d_all(filename[i], 'x') y = nc_read_from_file_2d_all(filename[i], 'y') # Save lai data as default for low and high vegetation lai_low = lai lai_high = lai # Read all tree parameters from file 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]) if (input_file_tree_crown_diameter[ii[i]] is not None): 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]) tree_crown_diameter = np.where( (tree_crown_diameter == 0.0) | (tree_crown_diameter == -1.0) ,fillvalues["tree_data"],tree_crown_diameter) else: tree_crown_diameter = np.ones((len(y),len(x))) tree_crown_diameter[:,:] = fillvalues["tree_data"] 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]) 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]) 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]) # Remove missing values from the data. Reasonable values will be set by the tree generator tree_height = np.where( (tree_height == 0.0) | (tree_height == -1.0) ,fillvalues["tree_data"],tree_height) tree_trunk_diameter = np.where( (tree_trunk_diameter == 0.0) | (tree_trunk_diameter == -1.0) ,fillvalues["tree_data"],tree_trunk_diameter) tree_type = np.where( (tree_type == 0.0) | (tree_type == -1.0) ,fillvalues["tree_data"],tree_type) patch_height = np.where( patch_height == -1.0 ,fillvalues["tree_data"],patch_height) # Convert trunk diameter from cm to m tree_trunk_diameter = np.where(tree_trunk_diameter != fillvalues["tree_data"], tree_trunk_diameter * 0.01,tree_trunk_diameter) # Temporarily change missing value for tree_type tree_type = np.where( (tree_type == fillvalues["tree_type"]),fillvalues["tree_data"],tree_type) # Compare patch height array with vegetation type and correct accordingly vegetation_type = nc_read_from_file_2d_all(filename[i], 'vegetation_type') # For zero-height patches, set vegetation_type to short grass and remove these pixels from the patch height array 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) 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) max_tree_height = max(tree_height.flatten()) max_patch_height = max(patch_height.flatten()) if ( (max_tree_height != fillvalues["tree_data"]) | (max_patch_height == fillvalues["tree_data"]) ): 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"]) # Remove LAD volumes that are inside buildings buildings_2d = nc_read_from_file_2d_all(filename[i], 'buildings_2d') for k in range(0,len(zlad)-1): lad[k,:,:] = np.where(buildings_2d == fillvalues["buildings_2d"],lad[k,:,:],fillvalues["tree_data"]) bad[k,:,:] = np.where(buildings_2d == fillvalues["buildings_2d"],bad[k,:,:],fillvalues["tree_data"]) tree_ids[k,:,:] = np.where(buildings_2d == fillvalues["buildings_2d"],tree_ids[k,:,:],fillvalues["tree_data"]) del buildings_2d nc_write_dimension(filename[i], 'zlad', zlad, datatypes["tree_data"]) nc_write_to_file_3d(filename[i], 'lad', lad, datatypes["tree_data"],'zlad','y','x',fillvalues["tree_data"]) nc_write_attribute(filename[i], 'lad', 'long_name', 'leaf area density') nc_write_attribute(filename[i], 'lad', 'units', '') nc_write_attribute(filename[i], 'lad', 'res_orig', domain_px[i]) nc_write_attribute(filename[i], 'lad', 'coordinates', 'E_UTM N_UTM lon lat') nc_write_attribute(filename[i], 'lad', 'grid_mapping', 'E_UTM N_UTM lon lat') nc_write_to_file_3d(filename[i], 'bad', bad, datatypes["tree_data"],'zlad','y','x',fillvalues["tree_data"]) nc_write_attribute(filename[i], 'bad', 'long_name', 'basal area density') nc_write_attribute(filename[i], 'bad', 'units', '') nc_write_attribute(filename[i], 'bad', 'res_orig', domain_px[i]) nc_write_attribute(filename[i], 'bad', 'coordinates', 'E_UTM N_UTM lon lat') nc_write_attribute(filename[i], 'bad', 'grid_mapping', 'E_UTM N_UTM lon lat') nc_write_to_file_3d(filename[i], 'tree_id', tree_ids, datatypes["tree_data"],'zlad','y','x',fillvalues["tree_data"]) nc_write_attribute(filename[i], 'tree_id', 'long_name', 'tree id') nc_write_attribute(filename[i], 'tree_id', 'units', '') nc_write_attribute(filename[i], 'tree_id', 'res_orig', domain_px[i]) nc_write_attribute(filename[i], 'tree_id', 'coordinates', 'E_UTM N_UTM lon lat') nc_write_attribute(filename[i], 'tree_id', 'grid_mapping', 'E_UTM N_UTM lon lat') del lai, lad, bad, tree_ids, zlad del vegetation_pars, tree_height, tree_crown_diameter, tree_trunk_diameter, tree_type, patch_height, x, y # Create vegetation patches for locations with high vegetation type for i in range(0,ndomains): if domain_canopy_patches[i]: # Load vegetation_type and lad array (at level z = 0) for re-processing vegetation_type = nc_read_from_file_2d_all(filename[i], 'vegetation_type') lad = nc_read_from_file_3d_all(filename[i], 'lad') zlad = nc_read_from_file_1d_all(filename[i], 'zlad') 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]) vegetation_pars = nc_read_from_file_3d_all(filename[i], 'vegetation_pars') lai = vegetation_pars[1,:,:] # 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 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"]) # Now, assign either the default LAI for high vegetation or the value stored in the LAI array. lai_high = np.where( (lai_high != fillvalues["tree_data"]) & (lai == fillvalues["tree_data"]), settings_lai_high_default, lai) # Define a patch height wherever it is missing, but where a high vegetation LAI was set patch_height = np.where ( (lai_high != fillvalues["tree_data"]) & (patch_height == fillvalues["tree_data"]), settings_patch_height_default, patch_height) # Remove pixels where street trees were already set patch_height = np.where ( (lad[0,:,:] != fillvalues["tree_data"]), fillvalues["tree_data"], patch_height) # Remove patch heights that have no lai_high value patch_height = np.where ( (lai_high == fillvalues["tree_data"]), fillvalues["tree_data"], patch_height) # For missing LAI values, set either the high vegetation default or the low vegetation default lai_high = np.where( (patch_height > 2.0) & (patch_height != fillvalues["tree_data"]) & (lai_high == fillvalues["tree_data"]),settings_lai_high_default,lai_high) lai_high = np.where( (patch_height <= 2.0) & (patch_height != fillvalues["tree_data"]) & (lai_high == fillvalues["tree_data"]),settings_lai_low_default,lai_high) if ( max(patch_height.flatten()) >= (2.0 * domain_dz[i]) ): print(" start calculating LAD (this might take some time)") lad_patch, patch_nz, status = process_patch(domain_dz[i],patch_height,max(zlad),lai_high,settings_lai_alpha,settings_lai_beta) 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,:,:] ) # Remove high vegetation wherever it is replaced by a leaf area density. This should effectively remove all high vegetation pixels vegetation_type = np.where((lad[0,:,:] != fillvalues["tree_data"]) & (vegetation_type != fillvalues["vegetation_type"]),3,vegetation_type) # If desired, remove all high vegetation. TODO: check if this is still necessary if not domain_high_vegetation[i]: 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) # Set default low LAI for pixels with an LAD (short grass below trees) lai_low = np.where( (lad[0,:,:] == fillvalues["tree_data"]), lai, settings_lai_low_default) # Fill low vegetation pixels without LAI set or with LAI = 0 with default value lai_low = np.where( ( (lai_low == fillvalues["tree_data"]) | (lai_low == 0.0) ) & (vegetation_type != fillvalues["vegetation_type"] ), settings_lai_low_default, lai_low) # Remove lai for pixels that have no vegetation_type lai_low = np.where( vegetation_type != fillvalues["vegetation_type"], lai_low, fillvalues["tree_data"]) # Overwrite lai in vegetation_parameters vegetation_pars[1,:,:] = np.copy(lai_low) nc_overwrite_to_file_3d(filename[i], 'vegetation_pars', vegetation_pars) # Overwrite lad array nc_overwrite_to_file_3d(filename[i], 'lad', lad) nc_overwrite_to_file_2d(filename[i], 'vegetation_type', vegetation_type) del vegetation_type, lad, lai, patch_height, vegetation_pars, zlad # Final consistency check for i in range(0,ndomains): vegetation_type = nc_read_from_file_2d_all(filename[i], 'vegetation_type') pavement_type = nc_read_from_file_2d_all(filename[i], 'pavement_type') building_type = nc_read_from_file_2d_all(filename[i], 'building_type') water_type = nc_read_from_file_2d_all(filename[i], 'water_type') soil_type = nc_read_from_file_2d_all(filename[i], '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"]) # Check for consistency and fill empty fields with default vegetation type consistency_array, test = check_consistency_3(vegetation_type,pavement_type,soil_type,fillvalues["vegetation_type"],fillvalues["pavement_type"],fillvalues["soil_type"]) surface_fraction = nc_read_from_file_3d_all(filename[i], 'surface_fraction') 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_overwrite_to_file_3d(filename[i], 'surface_fraction', surface_fraction) del vegetation_type, pavement_type, building_type, water_type, soil_type