#!/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 3567 2018-11-27 13:59:21Z sward $ # Initial revisions # # # # # # Description: # ------------ # Processing tool for creating PIDS conform static drivers from rastered NetCDF # input # # @Author Bjoern Maronga (maronga@muk.uni-hannover.de) # # @todo Remove high vegetation on demand # @todo Add vegetation_pars (LAI) # @todo Add building_pars (green roofs) # @todo Add LAD and BAD arrays (canopy generator) # @todo Make input files optional # @todo Allow for ASCII input of terrain height and building height # @todo Modularize reading config file #------------------------------------------------------------------------------# from palm_csd_files.palm_csd_netcdf_interface import * from palm_csd_files.palm_csd_tools import * import numpy as np def read_config_file(): import configparser from math import floor import numpy as np import os import subprocess as sub import sys # Check if configuration files exists and quit otherwise input_config = ".csd.config" for i in range(1,len(sys.argv)): input_config = str(sys.argv[i]) config = configparser.RawConfigParser(allow_no_value=True) if ( os.path.isfile(input_config) == False ): print ("Error. No configuration file " + input_config + " found.") raise SystemExit else: print(os.path.isfile(input_config)) config.read(input_config) # Definition of settings global settings_filename_out global settings_lai_season global settings_bridge_width global ndomains # Definition of global configuration parameters global global_acronym global global_angle global global_author global global_campaign global global_comment global global_contact global global_data_content global global_dependencies global global_institution global global_keywords global global_location global global_palm_version global global_references global global_site global global_source global global_version # Definition of domain parameters global domain_names global domain_px global domain_x0 global domain_y0 global domain_x1 global domain_y1 global domain_nx global domain_ny global domain_dz global domain_3d global domain_hv global domain_cg global domain_ip global domain_za global domain_parent # Definition of input data parameters global input_names global input_px global input_file_x global input_file_y global input_file_x_UTM global input_file_y_UTM global input_file_lat global input_file_lon global input_file_zt global input_file_buildings_2d global input_file_bridges_2d global input_file_building_id global input_file_bridges_id global input_file_building_type global input_file_building_type global input_file_vegetation_type global input_file_vegetation_height global input_file_pavement_type global input_file_water_type global input_file_street_type global input_file_street_crossings global input_file_soil_type global zt_all global zt_min settings_filename_out = "default" settings_lai_season = "summer" settings_bridge_width = 3.0 ndomains = 0 global_acronym = "" global_angle = "" global_author = "" global_campaign = "" global_comment = "" global_contact = "" global_data_content = "" global_dependencies = "" global_institution = "" global_keywords = "" global_location = "" global_palm_version = 6.0 global_references = "" global_site = "" global_source = "" global_version = 1 domain_names = [] domain_px = [] domain_x0 = [] domain_y0 = [] domain_x1 = [] domain_y1 = [] domain_nx = [] domain_ny = [] domain_dz = [] domain_3d = [] domain_hv = [] domain_cg = [] domain_ip = [] domain_za = [] domain_parent = [] zt_min = 0.0 zt_all = [] input_names = [] input_px = [] input_file_x = [] input_file_y = [] input_file_x_UTM = [] input_file_y_UTM = [] input_file_lat = [] input_file_lon = [] input_file_zt = [] input_file_buildings_2d = [] input_file_bridges_2d = [] input_file_building_id = [] input_file_bridges_id = [] input_file_building_type = [] input_file_vegetation_type = [] input_file_vegetation_height = [] input_file_pavement_type = [] input_file_water_type = [] input_file_street_type = [] input_file_street_crossings = [] input_file_soil_type = [] # Load all user parameters from config file for i in range(0,len(config.sections())): read_tmp = config.sections()[i] if ( read_tmp == 'global' ): global_acronym = config.get(read_tmp, 'acronym') global_angle = float(config.get(read_tmp, 'rotation_angle')) global_author = config.get(read_tmp, 'author') global_campaign = config.get(read_tmp, 'campaign') global_comment = config.get(read_tmp, 'comment') global_contact = config.get(read_tmp, 'contact_person') global_data_content = config.get(read_tmp, 'data_content') global_dependencies = config.get(read_tmp, 'dependencies') global_institution = config.get(read_tmp, 'institution') global_keywords = config.get(read_tmp, 'keywords') global_location = config.get(read_tmp, 'location') global_palm_version = float(config.get(read_tmp, 'palm_version')) global_references = config.get(read_tmp, 'references') global_site = config.get(read_tmp, 'site') global_source = config.get(read_tmp, 'source') global_version = float(config.get(read_tmp, 'version')) if ( read_tmp == 'settings' ): settings_filename_out = config.get(read_tmp, 'filename_out') settings_lai_season = config.get(read_tmp, 'lai_season') settings_bridge_width = float(config.get(read_tmp, 'bridge_width')) if ( read_tmp.split("_")[0] == 'domain' ): ndomains = ndomains + 1 domain_names.append(read_tmp.split("_")[1]) domain_px.append(float(config.get(read_tmp, 'pixel_size'))) domain_nx.append(int(config.get(read_tmp, 'nx'))) domain_ny.append(int(config.get(read_tmp, 'ny'))) domain_dz.append(float(config.get(read_tmp, 'dz'))) domain_3d.append(config.getboolean(read_tmp, 'buildings_3d')) domain_hv.append(config.getboolean(read_tmp, 'allow_high_vegetation')) domain_cg.append(config.getboolean(read_tmp, 'generate_vegetation_patches')) domain_ip.append(config.getboolean(read_tmp, 'interpolate_terrain')) domain_za.append(config.getboolean(read_tmp, 'use_palm_z_axis')) if domain_ip[ndomains-1] and not domain_za[ndomains-1]: domain_za[ndomains-1] = True print("+++ Overwrite user setting for use_palm_z_axis") domain_parent.append(config.get(read_tmp, 'domain_parent')) domain_x0.append(int(floor(float(config.get(read_tmp, 'origin_x'))/float(config.get(read_tmp, 'pixel_size'))))) domain_y0.append(int(floor(float(config.get(read_tmp, 'origin_y'))/float(config.get(read_tmp, 'pixel_size'))))) domain_x1.append(int(floor(float(config.get(read_tmp, 'origin_x'))/float(config.get(read_tmp, 'pixel_size'))) + int(config.get(read_tmp, 'nx')))) domain_y1.append(int(floor(float(config.get(read_tmp, 'origin_y'))/float(config.get(read_tmp, 'pixel_size'))) + int(config.get(read_tmp, 'ny')))) if ( read_tmp.split("_")[0] == 'input' ): input_names.append(read_tmp.split("_")[1]) input_px.append(float(config.get(read_tmp, 'pixel_size'))) input_file_x.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_x')) input_file_y.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_y')) input_file_lat.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_lat')) input_file_lon.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_lon')) input_file_x_UTM.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_x_UTM')) input_file_y_UTM.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_y_UTM')) input_file_zt.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_zt')) input_file_buildings_2d.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_buildings_2d')) input_file_bridges_2d.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_bridges_2d')) input_file_building_id.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_building_id')) input_file_bridges_id.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_bridges_id')) input_file_building_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_building_type')) input_file_vegetation_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_vegetation_type')) input_file_vegetation_height.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_vegetation_height')) input_file_pavement_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_pavement_type')) input_file_water_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_water_type')) input_file_street_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_street_type')) input_file_street_crossings.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_street_crossings')) #input_file_soil_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_soil_type')) return 0 ############################################################ # Start of main program datatypes = { "x": "f4", "y": "f4", "z": "f4", "lat": "f4", "lon": "f4", "E_UTM": "f4", "N_UTM": "f4", "zt": "f4", "buildings_2d": "f4", "buildings_3d": "b", "bridges_2d": "f4", "building_id": "i", "bridges_id": "i", "building_type": "b", "nsurface_fraction": "i", "vegetation_type": "b", "vegetation_height": "f4", "pavement_type": "b", "water_type": "b", "street_type": "b", "street_crossings": "b", "soil_type": "b", "surface_fraction": "f4" } fillvalues = { "lat": float(-9999.0), "lon": float(-9999.0), "E_UTM": float(-9999.0), "N_UTM": float(-9999.0), "zt": float(-9999.0), "buildings_2d": float(-9999.0), "buildings_3d": np.byte(-127), "bridges_2d": float(-9999.0), "building_id": int(-9999), "bridges_id": int(-9999), "building_type": np.byte(-127), "nsurface_fraction": int(-9999), "vegetation_type": np.byte(-127), "vegetation_height": float(-9999.0), "pavement_type": np.byte(-127), "water_type": np.byte(-127), "street_type": np.byte(-127), "street_crossings": np.byte(-127), "soil_type": np.byte(-127), "surface_fraction": float(-9999.0) } defaultvalues = { "lat": float(-9999.0), "lon": float(-9999.0), "E_UTM": float(-9999.0), "N_UTM": float(-9999.0), "zt": float(0.0), "buildings_2d": float(0.0), "buildings_3d": 0, "bridges_2d": float(0.0), "building_id": int(0), "bridges_id": int(0), "building_type": 1, "nsurface_fraction": int(-9999), "vegetation_type": 3, "vegetation_height": float(-9999.0), "pavement_type": 1, "water_type": 1, "street_type": 1, "street_crossings": 0, "soil_type": 1, "surface_fraction": float(0.0) } # Read configuration file and set parameters accordingly read_config_file() filename = [] ii = [] ii_parent = [] # Define indices and filenames for all domains and create netCDF files for i in range(0,ndomains): # Calculate indices and input files ii.append(input_px.index(domain_px[i])) filename.append(settings_filename_out + "_" + domain_names[i]) if domain_parent[i] is not None: ii_parent.append(domain_names.index(domain_parent[i])) else: ii_parent.append(None) x_UTM = nc_read_from_file_2d(input_file_x[ii[i]], "Band1", domain_x0[i], domain_x0[i], domain_y0[i], domain_y0[i]) y_UTM = nc_read_from_file_2d(input_file_y[ii[i]], "Band1", domain_x0[i], domain_x0[i], domain_y0[i], domain_y0[i]) lat = nc_read_from_file_2d(input_file_lat[ii[i]], "Band1", domain_x0[i], domain_x0[i], domain_y0[i], domain_y0[i]) lon = nc_read_from_file_2d(input_file_lon[ii[i]], "Band1", domain_x0[i], domain_x0[i], domain_y0[i], domain_y0[i]) # Create NetCDF output file and set global attributes nc_create_file(filename[i]) nc_write_global_attributes(filename[i],float(x_UTM[0,0]),float(y_UTM[0,0]),float(lat[0,0]),float(lon[0,0]),"",global_acronym,global_angle,global_author,global_campaign,global_comment,global_contact,global_data_content,global_dependencies,global_institution,global_keywords,global_location,global_palm_version,global_references,global_site,global_source,global_version) # Process terrain height for i in range(0,ndomains): # Read and write terrain height (zt) zt = nc_read_from_file_2d(input_file_zt[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) # Final step: add zt array to the global array zt_all.append(zt) del zt # Calculate the global (all domains) minimum of the terrain height. This value will be substracted for all # data sets zt_min = min(zt_all[0].flatten()) for i in range(0,ndomains): zt_min = min(zt_min,min(zt_all[i].flatten())) del zt_all[:] for i in range(0,ndomains): # Read and write terrain height (zt) zt = nc_read_from_file_2d(input_file_zt[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) x = nc_read_from_file_1d(input_file_x[ii[i]], "x", domain_x0[i], domain_x1[i]) y = nc_read_from_file_1d(input_file_y[ii[i]], "y", domain_y0[i], domain_y1[i]) print( "Shift terrain height by -" + str(zt_min)) zt = zt - zt_min # If necessary, interpolate parent domain terrain height on child domain grid and blend the two if domain_ip[i]: tmp_x0 = int( domain_x0[i] * domain_px[i] / domain_px[ii_parent[i]] ) - 1 tmp_y0 = int( domain_y0[i] * domain_px[i] / domain_px[ii_parent[i]] ) - 1 tmp_x1 = int( domain_x1[i] * domain_px[i] / domain_px[ii_parent[i]] ) + 1 tmp_y1 = int( domain_y1[i] * domain_px[i] / domain_px[ii_parent[i]] ) + 1 tmp_x = nc_read_from_file_1d(input_file_x[ii_parent[i]], "x", tmp_x0, tmp_x1) tmp_y = nc_read_from_file_1d(input_file_y[ii_parent[i]], "y", tmp_y0, tmp_y1) zt_parent = nc_read_from_file_2d(input_file_zt[ii_parent[i]], 'Band1', tmp_x0, tmp_x1, tmp_y0, tmp_y1) print( "Shift terrain height by -" + str(zt_min)) zt_parent = zt_parent - zt_min # Interpolate array and bring to PALM grid of child domain zt_ip = interpolate_2d(zt_parent,tmp_x,tmp_y,x,y) zt_ip = bring_to_palm_grid(zt_ip,x,y,domain_dz[ii_parent[i]]) # Shift the child terrain height according to the parent mean terrain height zt = zt - np.mean(zt) + np.mean(zt_ip) # Blend over the parent and child terrain height within a radius of 50 px zt = blend_array_2d(zt,zt_ip,50) # Final step: add zt array to the global array zt_all.append(zt) del zt # Read and shift x and y coordinates, shift terrain height according to its minimum value and write all data # to file for i in range(0,ndomains): # Read horizontal grid variables from zt file and write them to output file x = nc_read_from_file_1d(input_file_x[ii[i]], "x", domain_x0[i], domain_x1[i]) y = nc_read_from_file_1d(input_file_y[ii[i]], "y", domain_y0[i], domain_y1[i]) x = x - min(x.flatten()) + domain_px[i]/2.0 y = y - min(y.flatten()) + domain_px[i]/2.0 nc_write_dimension(filename[i], 'x', x, datatypes["x"]) nc_write_dimension(filename[i], 'y', y, datatypes["y"]) nc_write_attribute(filename[i], 'x', 'long_name', 'x') nc_write_attribute(filename[i], 'x', 'standard_name','projection_x_coordinate') nc_write_attribute(filename[i], 'x', 'units', 'm') nc_write_attribute(filename[i], 'y', 'long_name', 'x') nc_write_attribute(filename[i], 'y', 'standard_name', 'projection_y_coordinate') nc_write_attribute(filename[i], 'y', 'units', 'm') lat = nc_read_from_file_2d(input_file_lat[ii[i]], "Band1", domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) lon = nc_read_from_file_2d(input_file_lon[ii[i]], "Band1", domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) nc_write_to_file_2d(filename[i], 'lat', lat, datatypes["lat"],'y','x',fillvalues["lat"]) nc_write_to_file_2d(filename[i], 'lon', lon, datatypes["lon"],'y','x',fillvalues["lon"]) nc_write_attribute(filename[i], 'lat', 'long_name', 'latitude') nc_write_attribute(filename[i], 'lat', 'standard_name','latitude') nc_write_attribute(filename[i], 'lat', 'units', 'degrees_north') nc_write_attribute(filename[i], 'lon', 'long_name', 'longitude') nc_write_attribute(filename[i], 'lon', 'standard_name','longitude') nc_write_attribute(filename[i], 'lon', 'units', 'degrees_east') x_UTM = nc_read_from_file_2d(input_file_x_UTM[ii[i]], "Band1", domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) y_UTM = nc_read_from_file_2d(input_file_y_UTM[ii[i]], "Band1", domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) nc_write_to_file_2d(filename[i], 'E_UTM', x_UTM, datatypes["E_UTM"],'y','x',fillvalues["E_UTM"]) nc_write_to_file_2d(filename[i], 'N_UTM', y_UTM, datatypes["N_UTM"],'y','x',fillvalues["N_UTM"]) nc_write_attribute(filename[i], 'E_UTM', 'long_name', 'easting') nc_write_attribute(filename[i], 'E_UTM', 'standard_name','projection_x_coorindate') nc_write_attribute(filename[i], 'E_UTM', 'units', 'm') nc_write_attribute(filename[i], 'N_UTM', 'long_name', 'northing') nc_write_attribute(filename[i], 'N_UTM', 'standard_name','projection_y_coorindate') nc_write_attribute(filename[i], 'N_UTM', 'units', 'm') nc_write_crs(filename[i]) # If necessary, bring terrain height to PALM's vertical grid. This is either forced by the user or implicitly # by using interpolation for a child domain if domain_za[i]: zt_all[i] = bring_to_palm_grid(zt_all[i],x,y,domain_dz[i]) nc_write_to_file_2d(filename[i], 'zt', zt_all[i], datatypes["zt"],'y','x',fillvalues["zt"]) nc_write_attribute(filename[i], 'zt', 'long_name', 'orography') nc_write_attribute(filename[i], 'zt', 'units', 'm') nc_write_attribute(filename[i], 'zt', 'res_orig', domain_px[i]) nc_write_attribute(filename[i], 'zt', 'coordinates', 'E_UTM N_UTM lon lat') nc_write_attribute(filename[i], 'zt', 'grid_mapping', 'E_UTM N_UTM lon lat') del zt_all # Process building height, id, and type for i in range(0,ndomains): buildings_2d = nc_read_from_file_2d(input_file_buildings_2d[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) building_id = nc_read_from_file_2d(input_file_building_id[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) building_type = nc_read_from_file_2d(input_file_building_type[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) building_type[building_type == 255] = fillvalues["building_type"] building_type = np.where(building_type < 1,defaultvalues["building_type"],building_type) check = check_arrays_2(buildings_2d,building_id,fillvalues["buildings_2d"],fillvalues["building_id"]) if not check: buildings_2d = np.where(building_id != fillvalues["building_id"],buildings_2d,fillvalues["buildings_2d"]) building_id = np.where(buildings_2d == fillvalues["buildings_2d"],fillvalues["building_id"],building_id) print("Data check #1 " + str(check_arrays_2(buildings_2d,building_id,fillvalues["buildings_2d"],fillvalues["building_id"]))) check = check_arrays_2(buildings_2d,building_type,fillvalues["buildings_2d"],fillvalues["building_type"]) if not check: building_type = np.where(buildings_2d == fillvalues["buildings_2d"],fillvalues["building_type"],building_type) building_type = np.where((building_type == fillvalues["building_type"]) & (buildings_2d != fillvalues["buildings_2d"]),defaultvalues["building_type"],building_type) print("Data check #2 " + str(check_arrays_2(buildings_2d,building_type,fillvalues["buildings_2d"],fillvalues["building_type"]))) nc_write_to_file_2d(filename[i], 'buildings_2d', buildings_2d, datatypes["buildings_2d"],'y','x',fillvalues["buildings_2d"]) nc_write_attribute(filename[i], 'buildings_2d', 'long_name', 'buildings') nc_write_attribute(filename[i], 'buildings_2d', 'units', 'm') nc_write_attribute(filename[i], 'buildings_2d', 'res_orig', domain_px[i]) nc_write_attribute(filename[i], 'buildings_2d', 'lod', 1) nc_write_attribute(filename[i], 'buildings_2d', 'coordinates', 'E_UTM N_UTM lon lat') nc_write_attribute(filename[i], 'buildings_2d', 'grid_mapping', 'E_UTM N_UTM lon lat') nc_write_to_file_2d(filename[i], 'building_id', building_id, datatypes["building_id"],'y','x',fillvalues["building_id"]) nc_write_attribute(filename[i], 'building_id', 'long_name', 'building id') nc_write_attribute(filename[i], 'building_id', 'units', '') nc_write_attribute(filename[i], 'building_id', 'res_orig', domain_px[i]) nc_write_attribute(filename[i], 'building_id', 'coordinates', 'E_UTM N_UTM lon lat') nc_write_attribute(filename[i], 'building_id', 'grid_mapping', 'E_UTM N_UTM lon lat') nc_write_to_file_2d(filename[i], 'building_type', building_type, datatypes["building_type"],'y','x',fillvalues["building_type"]) nc_write_attribute(filename[i], 'building_type', 'long_name', 'building type') nc_write_attribute(filename[i], 'building_type', 'units', '') nc_write_attribute(filename[i], 'building_type', 'res_orig', domain_px[i]) nc_write_attribute(filename[i], 'building_type', 'coordinates', 'E_UTM N_UTM lon lat') nc_write_attribute(filename[i], 'building_type', 'grid_mapping', 'E_UTM N_UTM lon lat') del buildings_2d del building_id del building_type # Create 3d buildings if necessary. In that course, read bridge objects and add them to building layer for i in range(0,ndomains): if domain_3d[i]: x = nc_read_from_file_2d_all(filename[i], 'x') y = nc_read_from_file_2d_all(filename[i], 'y') buildings_2d = nc_read_from_file_2d_all(filename[i], 'buildings_2d') building_id = nc_read_from_file_2d_all(filename[i], 'building_id') bridges_2d = nc_read_from_file_2d(input_file_bridges_2d[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) bridges_id = nc_read_from_file_2d(input_file_bridges_id[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) bridges_2d = np.where(bridges_2d == 0.0,fillvalues["bridges_2d"],bridges_2d) building_id = np.where(bridges_2d == fillvalues["bridges_2d"],building_id,bridges_id) if np.any(buildings_2d != fillvalues["buildings_2d"]): buildings_3d, z = make_3d_from_2d(buildings_2d,x,y,domain_dz[i]) if np.any(bridges_2d != fillvalues["bridges_2d"]): buildings_3d = make_3d_from_bridges_2d(buildings_3d,bridges_2d,x,y,domain_dz[i],settings_bridge_width,fillvalues["bridges_2d"]) else: print("Skipping creation of 3D bridges (no bridges in domain)") nc_write_dimension(filename[i], 'z', z, datatypes["z"]) nc_write_attribute(filename[i], 'z', 'long_name', 'z') nc_write_attribute(filename[i], 'z', 'units', 'm') nc_write_to_file_3d(filename[i], 'buildings_3d', buildings_3d, datatypes["buildings_3d"],'z','y','x',fillvalues["buildings_3d"]) nc_write_attribute(filename[i], 'buildings_3d', 'long_name', 'buildings 3d') nc_write_attribute(filename[i], 'buildings_3d', 'units', '') nc_write_attribute(filename[i], 'buildings_3d', 'res_orig', domain_px[i]) nc_write_attribute(filename[i], 'buildings_3d', 'lod', 2) del buildings_3d else: print("Skipping creation of 3D buildings (no buildings in domain)") del bridges_2d, bridges_id, building_id # Read vegetation type, water_type, pavement_type, soil_type and make fields consistent for i in range(0,ndomains): building_type = nc_read_from_file_2d_all(filename[i], 'building_type') vegetation_type = nc_read_from_file_2d(input_file_vegetation_type[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) vegetation_type[vegetation_type == 255] = fillvalues["vegetation_type"] vegetation_type = np.where((vegetation_type < 1) & (vegetation_type != fillvalues["vegetation_type"]),defaultvalues["vegetation_type"],vegetation_type) pavement_type = nc_read_from_file_2d(input_file_pavement_type[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) pavement_type[pavement_type == 255] = fillvalues["pavement_type"] pavement_type = np.where((pavement_type < 1) & (pavement_type != fillvalues["pavement_type"]),defaultvalues["pavement_type"],pavement_type) water_type = nc_read_from_file_2d(input_file_water_type[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) water_type[water_type == 255] = fillvalues["water_type"] water_type = np.where((water_type < 1) & (water_type != fillvalues["water_type"]),defaultvalues["water_type"],water_type) # to do: replace by real soil input data soil_type = nc_read_from_file_2d(input_file_vegetation_type[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) soil_type[soil_type == 255] = fillvalues["soil_type"] soil_type = np.where((soil_type < 1) & (soil_type != fillvalues["soil_type"]),defaultvalues["soil_type"],soil_type) # Make arrays consistent # #1 Set vegetation type to missing for pixel where a pavement type is set vegetation_type = np.where((vegetation_type != fillvalues["vegetation_type"]) & (pavement_type != fillvalues["pavement_type"]),fillvalues["vegetation_type"],vegetation_type) # #2 Set vegetation type to missing for pixel where a building type is set vegetation_type = np.where((vegetation_type != fillvalues["vegetation_type"]) & (building_type != fillvalues["building_type"]) ,fillvalues["vegetation_type"],vegetation_type) # #3 Set vegetation type to missing for pixel where a building type is set vegetation_type = np.where((vegetation_type != fillvalues["vegetation_type"]) & (water_type != fillvalues["water_type"]),fillvalues["vegetation_type"],vegetation_type) # #4 Remove pavement for pixels with buildings pavement_type = np.where((pavement_type != fillvalues["pavement_type"]) & (building_type != fillvalues["building_type"]),fillvalues["pavement_type"],pavement_type) # #5 Remove pavement for pixels with water pavement_type = np.where((pavement_type != fillvalues["pavement_type"]) & (water_type != fillvalues["water_type"]),fillvalues["pavement_type"],pavement_type) # #6 Remove water for pixels with buildings water_type = np.where((water_type != fillvalues["water_type"]) & (building_type != fillvalues["building_type"]),fillvalues["water_type"],water_type) # #7 to be removed: set default soil type everywhere soil_type = np.where((vegetation_type != fillvalues["vegetation_type"]) | (pavement_type != fillvalues["pavement_type"]),defaultvalues["soil_type"],fillvalues["soil_type"]) # Check for consistency and fill empty fields with default vegetation type consistency_array, test = check_consistency_4(vegetation_type,building_type,pavement_type,water_type,fillvalues["vegetation_type"],fillvalues["building_type"],fillvalues["pavement_type"],fillvalues["water_type"]) if test: vegetation_type = np.where(consistency_array == 0,defaultvalues["vegetation_type"],vegetation_type) consistency_array, test = check_consistency_4(vegetation_type,building_type,pavement_type,water_type,fillvalues["vegetation_type"],fillvalues["building_type"],fillvalues["pavement_type"],fillvalues["water_type"]) # Create surface_fraction array x = nc_read_from_file_2d_all(filename[i], 'x') y = nc_read_from_file_2d_all(filename[i], 'y') nsurface_fraction = np.arange(0,3) surface_fraction = np.ones((len(nsurface_fraction),len(y),len(x))) surface_fraction[0,:,:] = np.where(vegetation_type != fillvalues["vegetation_type"], 1.0, 0.0) surface_fraction[1,:,:] = np.where(pavement_type != fillvalues["pavement_type"], 1.0, 0.0) surface_fraction[2,:,:] = np.where(water_type != fillvalues["water_type"], 1.0, 0.0) nc_write_dimension(filename[i], 'nsurface_fraction', nsurface_fraction, datatypes["nsurface_fraction"]) nc_write_to_file_3d(filename[i], 'surface_fraction', surface_fraction, datatypes["surface_fraction"],'nsurface_fraction','y','x',fillvalues["surface_fraction"]) nc_write_attribute(filename[i], 'surface_fraction', 'long_name', 'surface fraction') nc_write_attribute(filename[i], 'surface_fraction', 'units', '') nc_write_attribute(filename[i], 'surface_fraction', 'res_orig', domain_px[i]) del surface_fraction # Correct vegetation_type when a vegetation height is available and is indicative of low vegeetation vegetation_height = nc_read_from_file_2d(input_file_vegetation_height[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) vegetation_type = np.where((vegetation_height != fillvalues["vegetation_height"]) & (vegetation_height == 0.0) & ((vegetation_type == 4) | (vegetation_type == 7) | (vegetation_type == 17)), 3, vegetation_type) vegetation_height = np.where((vegetation_height != fillvalues["vegetation_height"]) & (vegetation_height == 0.0) & ((vegetation_type == 4) | (vegetation_type == 7) | (vegetation_type == 17)), fillvalues["vegetation_height"],vegetation_height) nc_write_to_file_2d(filename[i], 'vegetation_type', vegetation_type, datatypes["vegetation_type"],'y','x',fillvalues["vegetation_type"]) nc_write_attribute(filename[i], 'vegetation_type', 'long_name', 'vegetation type') nc_write_attribute(filename[i], 'vegetation_type', 'units', '') nc_write_attribute(filename[i], 'vegetation_type', 'res_orig', domain_px[i]) nc_write_attribute(filename[i], 'vegetation_type', 'coordinates', 'E_UTM N_UTM lon lat') nc_write_attribute(filename[i], 'vegetation_type', 'grid_mapping', 'E_UTM N_UTM lon lat') del vegetation_type nc_write_to_file_2d(filename[i], 'pavement_type', pavement_type, datatypes["pavement_type"],'y','x',fillvalues["pavement_type"]) nc_write_attribute(filename[i], 'pavement_type', 'long_name', 'pavement type') nc_write_attribute(filename[i], 'pavement_type', 'units', '') nc_write_attribute(filename[i], 'pavement_type', 'res_orig', domain_px[i]) nc_write_attribute(filename[i], 'pavement_type', 'coordinates', 'E_UTM N_UTM lon lat') nc_write_attribute(filename[i], 'pavement_type', 'grid_mapping', 'E_UTM N_UTM lon lat') del pavement_type nc_write_to_file_2d(filename[i], 'water_type', water_type, datatypes["water_type"],'y','x',fillvalues["water_type"]) nc_write_attribute(filename[i], 'water_type', 'long_name', 'water type') nc_write_attribute(filename[i], 'water_type', 'units', '') nc_write_attribute(filename[i], 'water_type', 'res_orig', domain_px[i]) nc_write_attribute(filename[i], 'water_type', 'coordinates', 'E_UTM N_UTM lon lat') nc_write_attribute(filename[i], 'water_type', 'grid_mapping', 'E_UTM N_UTM lon lat') del water_type nc_write_to_file_2d(filename[i], 'soil_type', soil_type, datatypes["soil_type"],'y','x',fillvalues["soil_type"]) nc_write_attribute(filename[i], 'soil_type', 'long_name', 'soil type') nc_write_attribute(filename[i], 'soil_type', 'units', '') nc_write_attribute(filename[i], 'soil_type', 'res_orig', domain_px[i]) nc_write_attribute(filename[i], 'soil_type', 'coordinates', 'E_UTM N_UTM lon lat') nc_write_attribute(filename[i], 'soil_type', 'grid_mapping', 'E_UTM N_UTM lon lat') del soil_type # pixels with bridges get building_type = 7 = bridge. This does not change the _type setting for the under-bridge # area if domain_3d[i]: if np.any(building_type != fillvalues["building_type"]): bridges_2d = nc_read_from_file_2d(input_file_bridges_2d[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) bridges_2d = np.where(bridges_2d == 0.0,fillvalues["bridges_2d"],bridges_2d) building_type = np.where(bridges_2d != fillvalues["bridges_2d"],7,building_type) nc_overwrite_to_file_2d(filename[i], 'building_type', building_type) del building_type del bridges_2d # Read/Write street type and street crossings for i in range(0,ndomains): street_type = nc_read_from_file_2d(input_file_street_type[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) street_type[street_type == 255] = fillvalues["street_type"] street_type = np.where((street_type < 1) & (street_type != fillvalues["street_type"]),defaultvalues["street_type"],street_type) nc_write_to_file_2d(filename[i], 'street_type', street_type, datatypes["street_type"],'y','x',fillvalues["street_type"]) nc_write_attribute(filename[i], 'street_type', 'long_name', 'street type') nc_write_attribute(filename[i], 'street_type', 'units', '') nc_write_attribute(filename[i], 'street_type', 'res_orig', domain_px[i]) nc_write_attribute(filename[i], 'street_type', 'coordinates', 'E_UTM N_UTM lon lat') nc_write_attribute(filename[i], 'street_type', 'grid_mapping', 'E_UTM N_UTM lon lat') del street_type street_crossings = nc_read_from_file_2d(input_file_street_crossings[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) street_crossings[street_crossings == 255] = fillvalues["street_crossings"] street_crossings = np.where((street_crossings < 1) & (street_crossings != fillvalues["street_crossings"]),defaultvalues["street_crossings"],street_crossings) nc_write_to_file_2d(filename[i], 'street_crossings', street_crossings, datatypes["street_crossings"],'y','x',fillvalues["street_crossings"]) nc_write_attribute(filename[i], 'street_crossings', 'long_name', 'street crossings') nc_write_attribute(filename[i], 'street_crossings', 'units', '') nc_write_attribute(filename[i], 'street_crossings', 'res_orig', domain_px[i]) nc_write_attribute(filename[i], 'street_crossings', 'coordinates', 'E_UTM N_UTM lon lat') nc_write_attribute(filename[i], 'street_crossings', 'grid_mapping', 'E_UTM N_UTM lon lat') del street_crossings