Home


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

source: palm/trunk/SCRIPTS/palm_csd @ 3585

Last change on this file since 3585 was 3567, checked in by maronga, 6 years ago

added first version of palm_csd (preprocessing tool for creating static drivers)
  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 35.3 KB
Line 
1#!/usr/bin/env python3
2# -*- coding: utf-8 -*-
3#--------------------------------------------------------------------------------#
4# This file is part of the PALM model system.
5#
6# PALM is free software: you can redistribute it and/or modify it under the terms
7# of the GNU General Public License as published by the Free Software Foundation,
8# either version 3 of the License, or (at your option) any later version.
9#
10# PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12# A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13#
14# You should have received a copy of the GNU General Public License along with
15# PALM. If not, see <http://www.gnu.org/licenses/>.
16#
17# Copyright 1997-2018  Leibniz Universitaet Hannover
18#--------------------------------------------------------------------------------#
19#
20# Current revisions:
21# -----------------
22#
23#
24# Former revisions:
25# -----------------
26# $Id: palm_csd 3567 2018-11-27 13:59:21Z forkel $
27# Initial revisions
28#
29#
30#
31#
32#
33# Description:
34# ------------
35# Processing tool for creating PIDS conform static drivers from rastered NetCDF
36# input
37#
38# @Author Bjoern Maronga (maronga@muk.uni-hannover.de)
39#
40# @todo Remove high vegetation on demand
41# @todo Add vegetation_pars (LAI)
42# @todo Add building_pars (green roofs)
43# @todo Add LAD and BAD arrays (canopy generator)
44# @todo Make input files optional
45# @todo Allow for ASCII input of terrain height and building height
46# @todo Modularize reading config file
47#------------------------------------------------------------------------------#
48
49from palm_csd_files.palm_csd_netcdf_interface import *
50from palm_csd_files.palm_csd_tools import *
51import numpy as np
52
53   
54def read_config_file():
55   
56   import configparser
57   from math import floor
58   import numpy as np
59   import os
60   import subprocess as sub
61   import sys
62
63#  Check if configuration files exists and quit otherwise
64   input_config = ".csd.config"
65   for i in range(1,len(sys.argv)): 
66      input_config = str(sys.argv[i])
67       
68   config = configparser.RawConfigParser(allow_no_value=True)
69
70   if ( os.path.isfile(input_config) == False ):
71      print ("Error. No configuration file " + input_config + " found.")
72      raise SystemExit   
73   else:
74      print(os.path.isfile(input_config))
75
76   config.read(input_config)
77
78
79#  Definition of settings
80   global settings_filename_out
81   global settings_lai_season
82   global settings_bridge_width
83   global ndomains
84
85#  Definition of global configuration parameters
86   global global_acronym
87   global global_angle
88   global global_author
89   global global_campaign
90   global global_comment
91   global global_contact
92   global global_data_content 
93   global global_dependencies
94   global global_institution
95   global global_keywords
96   global global_location
97   global global_palm_version
98   global global_references
99   global global_site
100   global global_source
101   global global_version
102
103
104#  Definition of domain parameters
105   global domain_names
106   global domain_px
107   global domain_x0
108   global domain_y0
109   global domain_x1
110   global domain_y1
111   global domain_nx
112   global domain_ny
113   global domain_dz
114   global domain_3d
115   global domain_hv
116   global domain_cg
117   global domain_ip
118   global domain_za
119   global domain_parent
120
121#  Definition of input data parameters
122   global input_names
123   global input_px
124 
125       
126   global input_file_x
127   global input_file_y
128   global input_file_x_UTM
129   global input_file_y_UTM   
130   global input_file_lat
131   global input_file_lon 
132   global input_file_zt
133   global input_file_buildings_2d
134   global input_file_bridges_2d
135   global input_file_building_id
136   global input_file_bridges_id   
137   global input_file_building_type
138   global input_file_building_type
139   global input_file_vegetation_type
140   global input_file_vegetation_height
141   global input_file_pavement_type
142   global input_file_water_type
143   global input_file_street_type
144   global input_file_street_crossings
145   global input_file_soil_type
146
147
148   global zt_all
149   global zt_min
150
151   settings_filename_out = "default"
152   settings_lai_season = "summer"
153   settings_bridge_width = 3.0
154   ndomains = 0
155
156   global_acronym = ""
157   global_angle = ""
158   global_author = ""
159   global_campaign = ""
160   global_comment = ""
161   global_contact = ""
162   global_data_content = ""
163   global_dependencies = ""
164   global_institution = ""
165   global_keywords = ""
166   global_location = ""
167   global_palm_version = 6.0
168   global_references = ""
169   global_site = ""
170   global_source = ""
171   global_version = 1
172
173   domain_names = []
174   domain_px = []
175   domain_x0 = []
176   domain_y0 = []
177   domain_x1 = []
178   domain_y1 = []
179   domain_nx = []
180   domain_ny = []
181   domain_dz = []
182   domain_3d = []
183   domain_hv = []
184   domain_cg = []
185   domain_ip = []
186   domain_za = []
187   domain_parent = []
188
189   zt_min = 0.0
190   zt_all = []
191
192   input_names = []
193   input_px = []
194
195   input_file_x = []
196   input_file_y = [] 
197   input_file_x_UTM = []
198   input_file_y_UTM = [] 
199   input_file_lat = []
200   input_file_lon = [] 
201   
202   input_file_zt = []
203   input_file_buildings_2d = []
204   input_file_bridges_2d = []
205   input_file_building_id = []
206   input_file_bridges_id = []   
207   input_file_building_type = []
208   input_file_vegetation_type = []
209   input_file_vegetation_height = []   
210   input_file_pavement_type = []
211   input_file_water_type = []
212   input_file_street_type = [] 
213   input_file_street_crossings = []     
214   input_file_soil_type = []
215
216   
217# Load all user parameters from config file
218   for i in range(0,len(config.sections())):
219
220      read_tmp = config.sections()[i]
221     
222      if ( read_tmp == 'global' ):
223         global_acronym = config.get(read_tmp, 'acronym')
224         global_angle = float(config.get(read_tmp, 'rotation_angle'))
225         global_author = config.get(read_tmp, 'author')
226         global_campaign = config.get(read_tmp, 'campaign')
227         global_comment = config.get(read_tmp, 'comment')
228         global_contact = config.get(read_tmp, 'contact_person')
229         global_data_content = config.get(read_tmp, 'data_content')   
230         global_dependencies = config.get(read_tmp, 'dependencies')
231         global_institution = config.get(read_tmp, 'institution')
232         global_keywords = config.get(read_tmp, 'keywords')
233         global_location = config.get(read_tmp, 'location')
234         global_palm_version = float(config.get(read_tmp, 'palm_version'))
235         global_references = config.get(read_tmp, 'references')
236         global_site = config.get(read_tmp, 'site')
237         global_source = config.get(read_tmp, 'source')
238         global_version = float(config.get(read_tmp, 'version'))
239
240      if ( read_tmp == 'settings' ):
241         settings_filename_out = config.get(read_tmp, 'filename_out')
242         settings_lai_season = config.get(read_tmp, 'lai_season') 
243         settings_bridge_width = float(config.get(read_tmp, 'bridge_width'))     
244 
245      if ( read_tmp.split("_")[0] == 'domain' ):
246         ndomains = ndomains + 1
247         domain_names.append(read_tmp.split("_")[1])
248         domain_px.append(float(config.get(read_tmp, 'pixel_size')))
249         domain_nx.append(int(config.get(read_tmp, 'nx')))
250         domain_ny.append(int(config.get(read_tmp, 'ny')))
251         domain_dz.append(float(config.get(read_tmp, 'dz')))       
252         domain_3d.append(config.getboolean(read_tmp, 'buildings_3d'))
253         domain_hv.append(config.getboolean(read_tmp, 'allow_high_vegetation'))
254         domain_cg.append(config.getboolean(read_tmp, 'generate_vegetation_patches'))       
255         domain_ip.append(config.getboolean(read_tmp, 'interpolate_terrain'))
256         domain_za.append(config.getboolean(read_tmp, 'use_palm_z_axis')) 
257         if domain_ip[ndomains-1] and not domain_za[ndomains-1]:
258            domain_za[ndomains-1] = True
259            print("+++ Overwrite user setting for use_palm_z_axis")
260         
261         domain_parent.append(config.get(read_tmp, 'domain_parent'))
262
263         domain_x0.append(int(floor(float(config.get(read_tmp, 'origin_x'))/float(config.get(read_tmp, 'pixel_size')))))
264         domain_y0.append(int(floor(float(config.get(read_tmp, 'origin_y'))/float(config.get(read_tmp, 'pixel_size')))))       
265         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'))))
266         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'))))
267
268      if ( read_tmp.split("_")[0] == 'input' ):
269         input_names.append(read_tmp.split("_")[1])
270         input_px.append(float(config.get(read_tmp, 'pixel_size')))
271         input_file_x.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_x'))
272         input_file_y.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_y')) 
273         input_file_lat.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_lat'))
274         input_file_lon.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_lon'))   
275         input_file_x_UTM.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_x_UTM'))
276         input_file_y_UTM.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_y_UTM')) 
277         input_file_zt.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_zt'))
278         input_file_buildings_2d.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_buildings_2d'))
279         input_file_bridges_2d.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_bridges_2d'))
280         input_file_building_id.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_building_id')) 
281         input_file_bridges_id.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_bridges_id')) 
282         input_file_building_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_building_type')) 
283         input_file_vegetation_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_vegetation_type')) 
284         input_file_vegetation_height.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_vegetation_height'))         
285         input_file_pavement_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_pavement_type'))   
286         input_file_water_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_water_type'))
287         input_file_street_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_street_type'))   
288         input_file_street_crossings.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_street_crossings')) 
289         #input_file_soil_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_soil_type'))   
290   return 0
291
292
293############################################################
294
295# Start of main program
296datatypes = {
297   "x": "f4", 
298   "y": "f4", 
299   "z": "f4", 
300   "lat": "f4", 
301   "lon": "f4",   
302   "E_UTM": "f4", 
303   "N_UTM": "f4", 
304   "zt": "f4",
305   "buildings_2d": "f4",
306   "buildings_3d": "b",
307   "bridges_2d": "f4",
308   "building_id": "i",
309   "bridges_id": "i",
310   "building_type": "b",
311   "nsurface_fraction": "i",
312   "vegetation_type": "b",
313   "vegetation_height": "f4",
314   "pavement_type": "b", 
315   "water_type": "b", 
316   "street_type": "b", 
317   "street_crossings": "b",     
318   "soil_type": "b",
319   "surface_fraction": "f4"
320   }
321
322fillvalues = {
323   "lat": float(-9999.0),
324   "lon": float(-9999.0),
325   "E_UTM": float(-9999.0),
326   "N_UTM": float(-9999.0),   
327   "zt": float(-9999.0),
328   "buildings_2d": float(-9999.0),
329   "buildings_3d": np.byte(-127),
330   "bridges_2d": float(-9999.0),
331   "building_id": int(-9999),
332   "bridges_id": int(-9999),
333   "building_type": np.byte(-127),
334   "nsurface_fraction": int(-9999),
335   "vegetation_type": np.byte(-127),
336   "vegetation_height": float(-9999.0),
337   "pavement_type": np.byte(-127),
338   "water_type": np.byte(-127),
339   "street_type": np.byte(-127), 
340   "street_crossings": np.byte(-127),   
341   "soil_type": np.byte(-127),
342   "surface_fraction": float(-9999.0)
343   }
344
345defaultvalues = {
346   "lat": float(-9999.0),
347   "lon": float(-9999.0),
348   "E_UTM": float(-9999.0),
349   "N_UTM": float(-9999.0), 
350   "zt": float(0.0),
351   "buildings_2d": float(0.0),
352   "buildings_3d": 0,   
353   "bridges_2d": float(0.0),
354   "building_id": int(0),
355   "bridges_id": int(0),
356   "building_type": 1,
357   "nsurface_fraction": int(-9999),
358   "vegetation_type": 3,
359   "vegetation_height": float(-9999.0),
360   "pavement_type": 1, 
361   "water_type": 1, 
362   "street_type": 1, 
363   "street_crossings": 0,   
364   "soil_type": 1,
365   "surface_fraction": float(0.0)
366   }
367 
368# Read configuration file and set parameters accordingly
369read_config_file()
370
371
372filename = []
373ii = []
374ii_parent = []
375# Define indices and filenames for all domains and create netCDF files
376for i in range(0,ndomains):
377
378#  Calculate indices and input files
379   ii.append(input_px.index(domain_px[i]))
380   filename.append(settings_filename_out + "_" + domain_names[i])
381   if domain_parent[i] is not None:
382      ii_parent.append(domain_names.index(domain_parent[i]))
383   else:
384      ii_parent.append(None)
385
386
387   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])
388   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])
389   lat = nc_read_from_file_2d(input_file_lat[ii[i]], "Band1", domain_x0[i], domain_x0[i], domain_y0[i], domain_y0[i])
390   lon = nc_read_from_file_2d(input_file_lon[ii[i]], "Band1", domain_x0[i], domain_x0[i], domain_y0[i], domain_y0[i]) 
391   
392#  Create NetCDF output file and set global attributes
393   nc_create_file(filename[i])
394   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)
395
396   
397# Process terrain height
398for i in range(0,ndomains):
399
400#  Read and write terrain height (zt)
401   zt = nc_read_from_file_2d(input_file_zt[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])
402
403#  Final step: add zt array to the global array   
404   zt_all.append(zt)
405   del zt
406 
407#  Calculate the global (all domains) minimum of the terrain height. This value will be substracted for all
408#  data sets
409zt_min = min(zt_all[0].flatten()) 
410for i in range(0,ndomains):
411   zt_min = min(zt_min,min(zt_all[i].flatten()))
412   
413del zt_all[:]
414
415for i in range(0,ndomains):
416
417#  Read and write terrain height (zt)
418   zt = nc_read_from_file_2d(input_file_zt[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])
419   x = nc_read_from_file_1d(input_file_x[ii[i]], "x", domain_x0[i], domain_x1[i])
420   y = nc_read_from_file_1d(input_file_y[ii[i]], "y", domain_y0[i], domain_y1[i])
421
422   print( "Shift terrain height by -" + str(zt_min))
423   zt = zt - zt_min
424
425#  If necessary, interpolate parent domain terrain height on child domain grid and blend the two
426   if domain_ip[i]:
427      tmp_x0 = int( domain_x0[i] * domain_px[i] / domain_px[ii_parent[i]] ) - 1
428      tmp_y0 = int( domain_y0[i] * domain_px[i] / domain_px[ii_parent[i]] ) - 1
429      tmp_x1 = int( domain_x1[i] * domain_px[i] / domain_px[ii_parent[i]] ) + 1
430      tmp_y1 = int( domain_y1[i] * domain_px[i] / domain_px[ii_parent[i]] ) + 1
431               
432      tmp_x = nc_read_from_file_1d(input_file_x[ii_parent[i]], "x", tmp_x0, tmp_x1)
433      tmp_y = nc_read_from_file_1d(input_file_y[ii_parent[i]], "y", tmp_y0, tmp_y1)
434
435      zt_parent = nc_read_from_file_2d(input_file_zt[ii_parent[i]], 'Band1', tmp_x0, tmp_x1, tmp_y0, tmp_y1)
436
437      print( "Shift terrain height by -" + str(zt_min))
438      zt_parent = zt_parent - zt_min
439
440#     Interpolate array and bring to PALM grid of child domain
441      zt_ip = interpolate_2d(zt_parent,tmp_x,tmp_y,x,y)
442      zt_ip = bring_to_palm_grid(zt_ip,x,y,domain_dz[ii_parent[i]])
443
444#     Shift the child terrain height according to the parent mean terrain height       
445      zt = zt - np.mean(zt)  + np.mean(zt_ip)
446 
447 
448#     Blend over the parent and child terrain height within a radius of 50 px 
449      zt = blend_array_2d(zt,zt_ip,50)
450   
451#  Final step: add zt array to the global array   
452   zt_all.append(zt)
453   del zt
454 
455
456# Read and shift x and y coordinates, shift terrain height according to its minimum value and write all data
457# to file 
458for i in range(0,ndomains):   
459#  Read horizontal grid variables from zt file and write them to output file
460   x = nc_read_from_file_1d(input_file_x[ii[i]], "x", domain_x0[i], domain_x1[i])
461   y = nc_read_from_file_1d(input_file_y[ii[i]], "y", domain_y0[i], domain_y1[i]) 
462   x = x - min(x.flatten()) + domain_px[i]/2.0
463   y = y - min(y.flatten()) + domain_px[i]/2.0 
464   nc_write_dimension(filename[i], 'x', x, datatypes["x"])
465   nc_write_dimension(filename[i], 'y', y, datatypes["y"]) 
466   nc_write_attribute(filename[i], 'x', 'long_name', 'x')
467   nc_write_attribute(filename[i], 'x', 'standard_name','projection_x_coordinate')   
468   nc_write_attribute(filename[i], 'x', 'units', 'm')
469   nc_write_attribute(filename[i], 'y', 'long_name', 'x')
470   nc_write_attribute(filename[i], 'y', 'standard_name', 'projection_y_coordinate')   
471   nc_write_attribute(filename[i], 'y', 'units', 'm')
472
473   lat = nc_read_from_file_2d(input_file_lat[ii[i]], "Band1", domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])
474   lon = nc_read_from_file_2d(input_file_lon[ii[i]], "Band1", domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) 
475
476   nc_write_to_file_2d(filename[i], 'lat', lat, datatypes["lat"],'y','x',fillvalues["lat"])
477   nc_write_to_file_2d(filename[i], 'lon', lon, datatypes["lon"],'y','x',fillvalues["lon"]) 
478   
479   nc_write_attribute(filename[i], 'lat', 'long_name', 'latitude')
480   nc_write_attribute(filename[i], 'lat', 'standard_name','latitude')   
481   nc_write_attribute(filename[i], 'lat', 'units', 'degrees_north')
482 
483   nc_write_attribute(filename[i], 'lon', 'long_name', 'longitude')
484   nc_write_attribute(filename[i], 'lon', 'standard_name','longitude')   
485   nc_write_attribute(filename[i], 'lon', 'units', 'degrees_east') 
486   
487   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])
488   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]) 
489   
490   nc_write_to_file_2d(filename[i], 'E_UTM', x_UTM, datatypes["E_UTM"],'y','x',fillvalues["E_UTM"])
491   nc_write_to_file_2d(filename[i], 'N_UTM', y_UTM, datatypes["N_UTM"],'y','x',fillvalues["N_UTM"]) 
492   
493   nc_write_attribute(filename[i], 'E_UTM', 'long_name', 'easting')
494   nc_write_attribute(filename[i], 'E_UTM', 'standard_name','projection_x_coorindate')   
495   nc_write_attribute(filename[i], 'E_UTM', 'units', 'm')
496 
497   nc_write_attribute(filename[i], 'N_UTM', 'long_name', 'northing')
498   nc_write_attribute(filename[i], 'N_UTM', 'standard_name','projection_y_coorindate')   
499   nc_write_attribute(filename[i], 'N_UTM', 'units', 'm') 
500 
501   nc_write_crs(filename[i])
502   
503 
504 
505#  If necessary, bring terrain height to PALM's vertical grid. This is either forced by the user or implicitly
506#  by using interpolation for a child domain       
507   if domain_za[i]:
508      zt_all[i] = bring_to_palm_grid(zt_all[i],x,y,domain_dz[i])
509
510   nc_write_to_file_2d(filename[i], 'zt', zt_all[i], datatypes["zt"],'y','x',fillvalues["zt"])
511   nc_write_attribute(filename[i], 'zt', 'long_name', 'orography')
512   nc_write_attribute(filename[i], 'zt', 'units', 'm')
513   nc_write_attribute(filename[i], 'zt', 'res_orig', domain_px[i]) 
514   nc_write_attribute(filename[i], 'zt', 'coordinates', 'E_UTM N_UTM lon lat') 
515   nc_write_attribute(filename[i], 'zt', 'grid_mapping', 'E_UTM N_UTM lon lat')   
516
517del zt_all
518
519
520#  Process building height, id, and type
521for i in range(0,ndomains): 
522   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])
523   
524   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]) 
525 
526   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]) 
527   building_type[building_type == 255] = fillvalues["building_type"]
528   building_type = np.where(building_type < 1,defaultvalues["building_type"],building_type)
529 
530   check = check_arrays_2(buildings_2d,building_id,fillvalues["buildings_2d"],fillvalues["building_id"])
531   if not check:
532      buildings_2d = np.where(building_id != fillvalues["building_id"],buildings_2d,fillvalues["buildings_2d"])
533      building_id = np.where(buildings_2d == fillvalues["buildings_2d"],fillvalues["building_id"],building_id)
534      print("Data check #1 " + str(check_arrays_2(buildings_2d,building_id,fillvalues["buildings_2d"],fillvalues["building_id"])))
535
536   check = check_arrays_2(buildings_2d,building_type,fillvalues["buildings_2d"],fillvalues["building_type"])
537   if not check:
538      building_type = np.where(buildings_2d == fillvalues["buildings_2d"],fillvalues["building_type"],building_type)
539      building_type = np.where((building_type == fillvalues["building_type"]) & (buildings_2d != fillvalues["buildings_2d"]),defaultvalues["building_type"],building_type)
540      print("Data check #2 " + str(check_arrays_2(buildings_2d,building_type,fillvalues["buildings_2d"],fillvalues["building_type"])))
541 
542   nc_write_to_file_2d(filename[i], 'buildings_2d', buildings_2d, datatypes["buildings_2d"],'y','x',fillvalues["buildings_2d"])
543   nc_write_attribute(filename[i], 'buildings_2d', 'long_name', 'buildings')
544   nc_write_attribute(filename[i], 'buildings_2d', 'units', 'm')
545   nc_write_attribute(filename[i], 'buildings_2d', 'res_orig', domain_px[i]) 
546   nc_write_attribute(filename[i], 'buildings_2d', 'lod', 1) 
547   nc_write_attribute(filename[i], 'buildings_2d', 'coordinates', 'E_UTM N_UTM lon lat') 
548   nc_write_attribute(filename[i], 'buildings_2d', 'grid_mapping', 'E_UTM N_UTM lon lat')   
549   
550   nc_write_to_file_2d(filename[i], 'building_id', building_id, datatypes["building_id"],'y','x',fillvalues["building_id"])
551   nc_write_attribute(filename[i], 'building_id', 'long_name', 'building id')
552   nc_write_attribute(filename[i], 'building_id', 'units', '')
553   nc_write_attribute(filename[i], 'building_id', 'res_orig', domain_px[i]) 
554   nc_write_attribute(filename[i], 'building_id', 'coordinates', 'E_UTM N_UTM lon lat') 
555   nc_write_attribute(filename[i], 'building_id', 'grid_mapping', 'E_UTM N_UTM lon lat')   
556   
557   nc_write_to_file_2d(filename[i], 'building_type', building_type, datatypes["building_type"],'y','x',fillvalues["building_type"])
558   nc_write_attribute(filename[i], 'building_type', 'long_name', 'building type')
559   nc_write_attribute(filename[i], 'building_type', 'units', '')
560   nc_write_attribute(filename[i], 'building_type', 'res_orig', domain_px[i]) 
561   nc_write_attribute(filename[i], 'building_type', 'coordinates', 'E_UTM N_UTM lon lat') 
562   nc_write_attribute(filename[i], 'building_type', 'grid_mapping', 'E_UTM N_UTM lon lat')   
563 
564del buildings_2d
565del building_id
566del building_type
567
568# Create 3d buildings if necessary. In that course, read bridge objects and add them to building layer
569for i in range(0,ndomains): 
570   
571   if domain_3d[i]:
572      x = nc_read_from_file_2d_all(filename[i], 'x')
573      y = nc_read_from_file_2d_all(filename[i], 'y') 
574      buildings_2d = nc_read_from_file_2d_all(filename[i], 'buildings_2d') 
575      building_id = nc_read_from_file_2d_all(filename[i], 'building_id') 
576     
577      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])
578      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])
579     
580      bridges_2d = np.where(bridges_2d == 0.0,fillvalues["bridges_2d"],bridges_2d)
581      building_id = np.where(bridges_2d == fillvalues["bridges_2d"],building_id,bridges_id)
582     
583     
584      if np.any(buildings_2d != fillvalues["buildings_2d"]):
585         buildings_3d, z = make_3d_from_2d(buildings_2d,x,y,domain_dz[i])
586         if np.any(bridges_2d != fillvalues["bridges_2d"]):
587            buildings_3d = make_3d_from_bridges_2d(buildings_3d,bridges_2d,x,y,domain_dz[i],settings_bridge_width,fillvalues["bridges_2d"])
588         else:
589            print("Skipping creation of 3D bridges (no bridges in domain)")
590           
591           
592         nc_write_dimension(filename[i], 'z', z, datatypes["z"]) 
593         nc_write_attribute(filename[i], 'z', 'long_name', 'z') 
594         nc_write_attribute(filename[i], 'z', 'units', 'm')
595 
596         nc_write_to_file_3d(filename[i], 'buildings_3d', buildings_3d, datatypes["buildings_3d"],'z','y','x',fillvalues["buildings_3d"])   
597         nc_write_attribute(filename[i], 'buildings_3d', 'long_name', 'buildings 3d')
598         nc_write_attribute(filename[i], 'buildings_3d', 'units', '')
599         nc_write_attribute(filename[i], 'buildings_3d', 'res_orig', domain_px[i]) 
600         nc_write_attribute(filename[i], 'buildings_3d', 'lod', 2) 
601         
602         del buildings_3d
603         
604      else:
605         print("Skipping creation of 3D buildings (no buildings in domain)")
606
607
608      del bridges_2d, bridges_id, building_id
609
610
611
612# Read vegetation type, water_type, pavement_type, soil_type and make fields consistent
613for i in range(0,ndomains): 
614
615   building_type = nc_read_from_file_2d_all(filename[i], 'building_type') 
616     
617   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])   
618   vegetation_type[vegetation_type == 255] = fillvalues["vegetation_type"]
619   vegetation_type = np.where((vegetation_type < 1) & (vegetation_type != fillvalues["vegetation_type"]),defaultvalues["vegetation_type"],vegetation_type)
620   
621   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])   
622   pavement_type[pavement_type == 255] = fillvalues["pavement_type"]
623   pavement_type = np.where((pavement_type < 1) & (pavement_type != fillvalues["pavement_type"]),defaultvalues["pavement_type"],pavement_type)
624
625   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])   
626   water_type[water_type == 255] = fillvalues["water_type"]
627   water_type = np.where((water_type < 1) & (water_type != fillvalues["water_type"]),defaultvalues["water_type"],water_type)
628 
629#  to do: replace by real soil input data
630   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]) 
631   soil_type[soil_type == 255] = fillvalues["soil_type"]
632   soil_type = np.where((soil_type < 1) & (soil_type != fillvalues["soil_type"]),defaultvalues["soil_type"],soil_type)
633
634#  Make arrays consistent
635#  #1 Set vegetation type to missing for pixel where a pavement type is set
636   vegetation_type = np.where((vegetation_type != fillvalues["vegetation_type"]) & (pavement_type != fillvalues["pavement_type"]),fillvalues["vegetation_type"],vegetation_type)
637
638#  #2 Set vegetation type to missing for pixel where a building type is set
639   vegetation_type = np.where((vegetation_type != fillvalues["vegetation_type"]) & (building_type != fillvalues["building_type"]) ,fillvalues["vegetation_type"],vegetation_type)
640
641#  #3 Set vegetation type to missing for pixel where a building type is set
642   vegetation_type = np.where((vegetation_type != fillvalues["vegetation_type"]) & (water_type != fillvalues["water_type"]),fillvalues["vegetation_type"],vegetation_type)   
643
644#  #4 Remove pavement for pixels with buildings
645   pavement_type = np.where((pavement_type != fillvalues["pavement_type"]) & (building_type != fillvalues["building_type"]),fillvalues["pavement_type"],pavement_type) 
646
647#  #5 Remove pavement for pixels with water
648   pavement_type = np.where((pavement_type != fillvalues["pavement_type"]) & (water_type != fillvalues["water_type"]),fillvalues["pavement_type"],pavement_type)   
649
650#  #6 Remove water for pixels with buildings
651   water_type = np.where((water_type != fillvalues["water_type"]) & (building_type != fillvalues["building_type"]),fillvalues["water_type"],water_type) 
652   
653
654#  #7 to be removed: set default soil type everywhere
655   soil_type = np.where((vegetation_type != fillvalues["vegetation_type"]) | (pavement_type != fillvalues["pavement_type"]),defaultvalues["soil_type"],fillvalues["soil_type"]) 
656 
657 
658#  Check for consistency and fill empty fields with default vegetation type
659   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"])
660   
661   if test:
662      vegetation_type = np.where(consistency_array == 0,defaultvalues["vegetation_type"],vegetation_type)
663      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"])
664   
665#  Create surface_fraction array
666   x = nc_read_from_file_2d_all(filename[i], 'x')
667   y = nc_read_from_file_2d_all(filename[i], 'y') 
668   nsurface_fraction = np.arange(0,3)
669   surface_fraction = np.ones((len(nsurface_fraction),len(y),len(x)))
670   
671   surface_fraction[0,:,:] = np.where(vegetation_type != fillvalues["vegetation_type"], 1.0, 0.0)
672   surface_fraction[1,:,:] = np.where(pavement_type != fillvalues["pavement_type"], 1.0, 0.0)
673   surface_fraction[2,:,:] = np.where(water_type != fillvalues["water_type"], 1.0, 0.0)
674   
675   nc_write_dimension(filename[i], 'nsurface_fraction', nsurface_fraction, datatypes["nsurface_fraction"]) 
676   nc_write_to_file_3d(filename[i], 'surface_fraction', surface_fraction, datatypes["surface_fraction"],'nsurface_fraction','y','x',fillvalues["surface_fraction"])     
677   nc_write_attribute(filename[i], 'surface_fraction', 'long_name', 'surface fraction')
678   nc_write_attribute(filename[i], 'surface_fraction', 'units', '')
679   nc_write_attribute(filename[i], 'surface_fraction', 'res_orig', domain_px[i])     
680   del surface_fraction   
681   
682 
683#  Correct vegetation_type when a vegetation height is available and is indicative of low vegeetation
684   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])   
685         
686   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)
687   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)
688
689
690   nc_write_to_file_2d(filename[i], 'vegetation_type', vegetation_type, datatypes["vegetation_type"],'y','x',fillvalues["vegetation_type"])     
691   nc_write_attribute(filename[i], 'vegetation_type', 'long_name', 'vegetation type')
692   nc_write_attribute(filename[i], 'vegetation_type', 'units', '')
693   nc_write_attribute(filename[i], 'vegetation_type', 'res_orig', domain_px[i]) 
694   nc_write_attribute(filename[i], 'vegetation_type', 'coordinates', 'E_UTM N_UTM lon lat') 
695   nc_write_attribute(filename[i], 'vegetation_type', 'grid_mapping', 'E_UTM N_UTM lon lat')   
696   del vegetation_type
697
698   nc_write_to_file_2d(filename[i], 'pavement_type', pavement_type, datatypes["pavement_type"],'y','x',fillvalues["pavement_type"])   
699   nc_write_attribute(filename[i], 'pavement_type', 'long_name', 'pavement type')
700   nc_write_attribute(filename[i], 'pavement_type', 'units', '')
701   nc_write_attribute(filename[i], 'pavement_type', 'res_orig', domain_px[i]) 
702   nc_write_attribute(filename[i], 'pavement_type', 'coordinates', 'E_UTM N_UTM lon lat') 
703   nc_write_attribute(filename[i], 'pavement_type', 'grid_mapping', 'E_UTM N_UTM lon lat')   
704   del pavement_type
705
706   nc_write_to_file_2d(filename[i], 'water_type', water_type, datatypes["water_type"],'y','x',fillvalues["water_type"]) 
707   nc_write_attribute(filename[i], 'water_type', 'long_name', 'water type')
708   nc_write_attribute(filename[i], 'water_type', 'units', '')
709   nc_write_attribute(filename[i], 'water_type', 'res_orig', domain_px[i]) 
710   nc_write_attribute(filename[i], 'water_type', 'coordinates', 'E_UTM N_UTM lon lat') 
711   nc_write_attribute(filename[i], 'water_type', 'grid_mapping', 'E_UTM N_UTM lon lat')   
712   del water_type
713
714   nc_write_to_file_2d(filename[i], 'soil_type', soil_type, datatypes["soil_type"],'y','x',fillvalues["soil_type"]) 
715   nc_write_attribute(filename[i], 'soil_type', 'long_name', 'soil type')
716   nc_write_attribute(filename[i], 'soil_type', 'units', '')
717   nc_write_attribute(filename[i], 'soil_type', 'res_orig', domain_px[i])   
718   nc_write_attribute(filename[i], 'soil_type', 'coordinates', 'E_UTM N_UTM lon lat') 
719   nc_write_attribute(filename[i], 'soil_type', 'grid_mapping', 'E_UTM N_UTM lon lat')   
720   del soil_type
721
722
723
724
725
726#  pixels with bridges get building_type = 7 = bridge. This does not change the _type setting for the under-bridge
727#  area
728   if domain_3d[i]:
729      if np.any(building_type != fillvalues["building_type"]):
730 
731         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])         
732         bridges_2d = np.where(bridges_2d == 0.0,fillvalues["bridges_2d"],bridges_2d)
733         building_type = np.where(bridges_2d != fillvalues["bridges_2d"],7,building_type)
734         nc_overwrite_to_file_2d(filename[i], 'building_type', building_type) 
735 
736         del building_type 
737         del bridges_2d
738
739# Read/Write street type and street crossings
740for i in range(0,ndomains): 
741
742   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])   
743   street_type[street_type == 255] = fillvalues["street_type"]
744   street_type = np.where((street_type < 1) & (street_type != fillvalues["street_type"]),defaultvalues["street_type"],street_type)
745   
746   nc_write_to_file_2d(filename[i], 'street_type', street_type, datatypes["street_type"],'y','x',fillvalues["street_type"]) 
747   nc_write_attribute(filename[i], 'street_type', 'long_name', 'street type')
748   nc_write_attribute(filename[i], 'street_type', 'units', '')
749   nc_write_attribute(filename[i], 'street_type', 'res_orig', domain_px[i]) 
750   nc_write_attribute(filename[i], 'street_type', 'coordinates', 'E_UTM N_UTM lon lat') 
751   nc_write_attribute(filename[i], 'street_type', 'grid_mapping', 'E_UTM N_UTM lon lat')   
752   del street_type
753   
754   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])   
755   street_crossings[street_crossings == 255] = fillvalues["street_crossings"]
756   street_crossings = np.where((street_crossings < 1) & (street_crossings != fillvalues["street_crossings"]),defaultvalues["street_crossings"],street_crossings)
757   
758   nc_write_to_file_2d(filename[i], 'street_crossings', street_crossings, datatypes["street_crossings"],'y','x',fillvalues["street_crossings"]) 
759   nc_write_attribute(filename[i], 'street_crossings', 'long_name', 'street crossings')
760   nc_write_attribute(filename[i], 'street_crossings', 'units', '')
761   nc_write_attribute(filename[i], 'street_crossings', 'res_orig', domain_px[i]) 
762   nc_write_attribute(filename[i], 'street_crossings', 'coordinates', 'E_UTM N_UTM lon lat') 
763   nc_write_attribute(filename[i], 'street_crossings', 'grid_mapping', 'E_UTM N_UTM lon lat')   
764   del street_crossings
Note: See TracBrowser for help on using the repository browser.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |