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 kanani $ |
---|
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 | |
---|
49 | from palm_csd_files.palm_csd_netcdf_interface import * |
---|
50 | from palm_csd_files.palm_csd_tools import * |
---|
51 | import numpy as np |
---|
52 | |
---|
53 | |
---|
54 | def 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 |
---|
296 | datatypes = { |
---|
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 | |
---|
322 | fillvalues = { |
---|
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 | |
---|
345 | defaultvalues = { |
---|
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 |
---|
369 | read_config_file() |
---|
370 | |
---|
371 | |
---|
372 | filename = [] |
---|
373 | ii = [] |
---|
374 | ii_parent = [] |
---|
375 | # Define indices and filenames for all domains and create netCDF files |
---|
376 | for 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 |
---|
398 | for 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 |
---|
409 | zt_min = min(zt_all[0].flatten()) |
---|
410 | for i in range(0,ndomains): |
---|
411 | zt_min = min(zt_min,min(zt_all[i].flatten())) |
---|
412 | |
---|
413 | del zt_all[:] |
---|
414 | |
---|
415 | for 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 |
---|
458 | for 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 | |
---|
517 | del zt_all |
---|
518 | |
---|
519 | |
---|
520 | # Process building height, id, and type |
---|
521 | for 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 | |
---|
564 | del buildings_2d |
---|
565 | del building_id |
---|
566 | del building_type |
---|
567 | |
---|
568 | # Create 3d buildings if necessary. In that course, read bridge objects and add them to building layer |
---|
569 | for 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 |
---|
613 | for 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 |
---|
740 | for 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 |
---|