source: palm/trunk/SCRIPTS/create_basic_static_driver.py @ 4841

Last change on this file since 4841 was 4796, checked in by gronemeier, 4 years ago

add an example to create_basic_static_driver.py which creates static driver of urban_environment test case; update static driver of urban_environment test case

  • Property svn:executable set to *
File size: 19.3 KB
RevLine 
[4796]1#!/usr/bin/env python3
2# -------------------------------------------------------------------- #
[4072]3# This file is part of the PALM model system.
4#
[4796]5# PALM is free software: you can redistribute it and/or modify it under
6# the terms of the GNU General Public License as published by the Free
7# Software Foundation, either version 3 of the License, or (at your
8# option) any later version.
[4072]9#
[4796]10# PALM is distributed in the hope that it will be useful, but WITHOUT
11# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
13# for more details.
[4072]14#
[4796]15# You should have received a copy of the GNU General Public License
16# along with PALM. If not, see <http://www.gnu.org/licenses/>.
[4072]17#
[4370]18# Copyright 1997-2020 Leibniz Universitaet Hannover
[4796]19# -------------------------------------------------------------------- #
[4072]20
[4796]21import datetime
[4072]22import math
[4796]23import pytz
24
25from netCDF4 import Dataset
[4072]26import numpy as np
27
28
29class StaticDriver:
[4796]30    """This is an example script to generate static drivers for PALM.
31
[4072]32    You can use it as a starting point for creating your setup specific
[4796]33    driver.
[4072]34    """
[4796]35
[4072]36    def __init__(self):
[4796]37        """Open the static driver as netCDF4 file. Here, you have to
38        give the full path to the static driver that shall be created.
39        Existing file with same name is deleted.
[4072]40        """
41        print('Opening file...')
[4796]42        self.nc_file = Dataset('example_static_file.nc', 'w', format='NETCDF4')
[4072]43
44    def write_global_attributes(self):
[4796]45        """Write global attributes to static driver."""
[4072]46        print("Writing global attributes...")
[4796]47
48        # Optional global attributes
49        # --------------------------
50        self.nc_file.title = 'Example PALM static driver'
51        self.nc_file.author = 'PALM user'
[4072]52        self.nc_file.institution = 'Institut of Meteorology and Climatology,' \
[4796]53            'Leibniz University Hannover'
54        self.nc_file.comment = 'Generic crossing example'
55        self.nc_file.creation_date = \
56            pytz.utc.localize(datetime.datetime.utcnow()).strftime('%y-%m-%d %H:%M:%S %z')[:-2]
[4072]57        self.nc_file.history = ''
[4796]58        self.nc_file.keywords = 'example, PALM-4U'
59        self.nc_file.license = ''
60        self.nc_file.palm_version = ''
61        self.nc_file.references = ''
62        self.nc_file.source = 'PALM trunk'
63        self.nc_file.version = '1'
[4072]64
[4796]65        # Mandatory global attributes
66        # ---------------------------
67        self.nc_file.Conventions = 'CF-1.7'
68        self.nc_file.origin_lat = 52.50965  # (overwrite initialization_parameters)
69        self.nc_file.origin_lon = 13.3139  # Used to initialize Coriolis parameter
70        self.nc_file.origin_time = '2019-03-06 10:00:00 +00'
71        self.nc_file.origin_x = 3455249.0
72        self.nc_file.origin_y = 5424815.0
73        self.nc_file.origin_z = 0.0
74        self.nc_file.rotation_angle = 0.0
[4072]75
76    def define_dimensions(self):
[4796]77        """Set dimensions on which variables are defined."""
[4072]78        print("Writing dimensions...")
[4796]79
80        # Specify general grid parameters
81        # These values must equal to those set in the initialization_parameters
82        self.nx = 19
83        self.ny = 19
84        self.nz = 20
85        dx = 2
86        dy = 2
87        dz = 2
88
89        # Create soil grid (only relevant if land surface module is used)
[4072]90        dz_soil = np.array((0.01, 0.02, 0.04, 0.06, 0.14, 0.26, 0.54, 1.86))
91        zsoil_fullLayers = np.zeros_like(dz_soil)
[4796]92        zsoil_fullLayers = np.around(
93            [np.sum(dz_soil[:zs]) for zs in np.arange(1, len(dz_soil)+1)], 2)
[4072]94        zsoil_array = zsoil_fullLayers - dz_soil/2.
[4796]95
96        # Coordinates
97        # -----------
98        self.nc_file.createDimension('x', self.nx+1)
99        self.x = self.nc_file.createVariable('x', 'f4', ('x',))
100        self.x.long_name = 'distance to origin in x-direction'
[4072]101        self.x.units = 'm'
[4796]102        self.x.axis = 'X'
103        self.x[:] = np.arange(0, (self.nx+1)*dx, dx) + 0.5 * dx
104
[4072]105        self.nc_file.createDimension('y', self.ny+1)
[4796]106        self.y = self.nc_file.createVariable('y', 'f4', ('y',))
107        self.y.long_name = 'distance to origin in y-direction'
[4072]108        self.y.units = 'm'
[4796]109        self.y.axis = 'Y'
110        self.y[:] = np.arange(0, (self.ny+1)*dy, dy) + 0.5 * dy
111
112        # NOTE if your simulation uses a stretched vertical grid, you
113        # need to modify the z coordinates, e.g. z_array = (...)
114        z_array = np.append(0, np.arange(dz/2, (self.nz)*dz, dz))
115        self.nc_file.createDimension('z', self.nz+1)
116        self.z = self.nc_file.createVariable('z', 'f4', ('z',))
117        self.z.long_name = 'height above origin'
[4072]118        self.z.units = 'm'
[4796]119        self.z.axis = 'Z'
120        self.z.positive = 'up'
[4072]121        self.z[:] = z_array
122
[4796]123        zlad_array = self.z[:6]
124        self.nc_file.createDimension('zlad', len(zlad_array))
125        self.zlad = self.nc_file.createVariable('zlad', 'f4', ('zlad',))
126        self.zlad.long_name = 'height above ground'
127        self.zlad.units = 'm'
128        self.zlad.axis = 'Z'
129        self.zlad.positive = 'up'
130        self.zlad[:] = zlad_array
[4072]131
[4796]132        # self.nc_file.createDimension('zsoil', len(zsoil_array))
133        # self.zsoil = self.nc_file.createVariable('zsoil', 'f4', ('zsoil',))
134        # self.zsoil.long_name = 'depth in the soil'
135        # self.zsoil.units = 'm'
136        # self.zsoil.axis = 'Z'
137        # self.zsoil.positive = 'down'
138        # self.zsoil[:] = zsoil_array
[4072]139
[4796]140        # Other dimensions
141        # ----------------
142        # self.nc_file.createDimension('nalbedo_pars', 3)
143        # self.nalbedo_pars = self.nc_file.createVariable(
144        #     'nalbedo_pars', 'i4', ('nalbedo_pars',))
145        # self.nalbedo_pars[:] = np.arange(8)
146        #
147        # self.nc_file.createDimension('nbuilding_pars', 46)
148        # self.nbuilding_pars = self.nc_file.createVariable(
149        #     'nbuilding_pars', 'i4', ('nbuilding_pars',))
150        # self.nbuilding_pars[:] = np.arange(46)
151        #
152        # self.nc_file.createDimension('npavement_pars', 4)
153        # self.npavement_pars = self.nc_file.createVariable(
154        #     'npavement_pars', 'i4', ('npavement_pars',))
155        # self.npavement_pars[:] = np.arange(4)
156        #
157        # self.nc_file.createDimension('npavement_subsurface_pars', 2)
158        # self.npavement_subsurface_pars = self.nc_file.createVariable(
159        #     'npavement_subsurface_pars', 'i4', ('npavement_subsurface_pars',))
160        # self.npavement_subsurface_pars[:] = np.arange(2)
161        #
162        # self.nc_file.createDimension('nsoil_pars', 8)
163        # self.nsoil_pars = self.nc_file.createVariable(
164        #     'nsoil_pars', 'i4', ('nsoil_pars',))
165        # self.nsoil_pars[:] = np.arange(8)
166        #
167        self.nc_file.createDimension('nsurface_fraction', 3)
168        self.nsurface_fraction = self.nc_file.createVariable(
169            'nsurface_fraction', 'i4', ('nsurface_fraction',))
170        self.nsurface_fraction[:] = np.arange(3)
171
172        # self.nc_file.createDimension('nvegetation_pars', 12)
173        # self.nvegetation_pars = self.nc_file.createVariable(
174        #     'nvegetation_pars', 'i4', ('nvegetation_pars',))
175        # self.nvegetation_pars[:] = np.arange(12)
176        #
177        # self.nc_file.createDimension('nwater_pars', 6)
178        # self.nwater_pars = self.nc_file.createVariable(
179        #     'nwater_pars', 'i4', ('nwater_pars',))
180        # self.nwater_pars[:] = np.arange(7)
181
182    def define_variables(self):
183        """Define variables for the static driver.
184
185        Be aware that some variables depend on others. For a description
186        of each variable, please have a look at the documentation at
[4072]187        palm.muk.uni-hannover.de/trac/wiki/doc/app/iofiles/pids/static
[4796]188
[4072]189        An example of how you modify the variables is given below:
[4796]190
191        building_2d_array = np.ones((self.ny+1, self.nx+1)) * -9999.0
[4072]192        south_wall, north_wall, left_wall, right_wall = 20, 25, 20, 25
[4796]193        building_2d_array[
194            south_wall:north_wall,
195            left_wall:right_wall
196            ] = 50
[4072]197        nc_buildings_2d = self.nc_file.createVariable(
[4796]198            'buildings_2d', 'f4', ('y', 'x'), fill_value=-9999.0)
[4072]199        nc_buildings_2d.lod = 1
[4796]200        nc_buildings_2d[:, :] = building_2d_array
[4072]201        """
202        print("Writing variables...")
[4796]203
204        # Topography set-up
205        # -----------------
206        nc_building_id = self.nc_file.createVariable(
207            'building_id', 'i4', ('y', 'x'), fill_value=-9999)
208        nc_building_id.long_name = "building id number"
209        nc_building_id.units = "1"
210        nc_building_id[:, :] = nc_building_id._FillValue
211
212        nc_buildings_2d = self.nc_file.createVariable(
213            'buildings_2d', 'f4', ('y', 'x'), fill_value=-9999.0)
214        nc_buildings_2d.long_name = "building height"
215        nc_buildings_2d.units = "m"
216        nc_buildings_2d.lod = np.int32(1)
217        nc_buildings_2d[:, :] = nc_buildings_2d._FillValue
218
219        nc_buildings_3d = self.nc_file.createVariable(
220            'buildings_3d', 'i1', ('z', 'y', 'x'), fill_value=-127)
221        nc_buildings_3d.long_name = "building flag"
222        nc_buildings_3d.units = "1"
223        nc_buildings_3d.lod = np.int32(2)
224        nc_buildings_3d[:, :, :] = nc_buildings_3d._FillValue
225
226        nc_zt = self.nc_file.createVariable(
227            'zt', 'f4', ('y', 'x'), fill_value=-9999.0)
228        nc_zt.long_name = 'terrain height'
229        nc_zt.units = 'm'
230        nc_zt[:, :] = nc_zt._FillValue
231
232        # nc_z0 = self.nc_file.createVariable(
233        #     'z0', 'f4', ('y', 'x'), fill_value=-9999.0)
234        # nc_z0.long_name = 'roughness length for momentum'
235        # nc_z0.units = 'm'
236        # nc_z0[:, :] = nc_z0._FillValue
237
238        # Main surface clasification
239        # --------------------------
240        # nc_albedo_type = self.nc_file.createVariable(
241        #     'albedo_type', 'f4', ('y', 'x'), fill_value=-9999.0)
242        # nc_albedo_type.long_name = "albedo type"
243        # nc_albedo_type.units = "1"
244        # nc_albedo_type[:, :] = nc_albedo_type._FillValue
245        #
246        nc_building_type = self.nc_file.createVariable(
247            'building_type', 'i1', ('y', 'x'), fill_value=-127)
248        nc_building_type.long_name = "building type classification"
249        nc_building_type.units = "1"
250        nc_building_type[:, :] = nc_building_type._FillValue
251
252        nc_pavement_type = self.nc_file.createVariable(
253            'pavement_type', 'i1', ('y', 'x'), fill_value=-127)
254        nc_pavement_type.long_name = "pavement type classification"
255        nc_pavement_type.units = "1"
256        nc_pavement_type[:, :] = nc_pavement_type._FillValue
257
258        nc_soil_type = self.nc_file.createVariable(
259            'soil_type', 'i1', ('y', 'x'), fill_value=-127)
260        nc_soil_type.long_name = "soil type classification"
261        nc_soil_type.units = "1"
262        nc_soil_type.lod = np.int32(1)
263        nc_soil_type[:, :] = nc_soil_type._FillValue
264
265        # NOTE alternatively, define different soil types for different
266        # soil layers
267        # nc_soil_type = self.nc_file.createVariable(
268        #     'soil_type', 'i1', ('zsoil', 'y', 'x'), fill_value=-127)
269        # nc_soil_type.long_name = "soil type"
270        # nc_soil_type.units = "1"
271        # nc_soil_type.lod = np.int32(2)
272        # nc_soil_type[:, :, :] = nc_soil_type._FillValue
273        #
274        # nc_street_crossing = self.nc_file.createVariable(
275        #     'street_crossing', 'i1', ('y', 'x'),
276        #     fill_value=-127)
277        # nc_street_crossing.long_name = "street crossing"
278        # nc_street_crossing.units = "1"
279        # nc_street_crossing.valid_range = np.byte(1)
280        # nc_street_crossing[:, :] = nc_street_crossing._FillValue
281        #
282        nc_street_type = self.nc_file.createVariable(
283            'street_type', 'i1', ('y', 'x'),
284            fill_value=-127)
285        nc_street_type.long_name = "street type classification"
286        nc_street_type.units = "1"
287        nc_street_type[:, :] = nc_street_type._FillValue
288
289        nc_surface_fraction = self.nc_file.createVariable(
290            'surface_fraction', 'f4', ('nsurface_fraction', 'y', 'x'), fill_value=-9999.0)
291        nc_surface_fraction.long_name = "surface fraction"
292        nc_surface_fraction.units = "1"
293        nc_surface_fraction[0, :, :] = nc_surface_fraction._FillValue  # vegetation fraction
294        nc_surface_fraction[1, :, :] = nc_surface_fraction._FillValue  # pavement fraction
295        nc_surface_fraction[2, :, :] = nc_surface_fraction._FillValue  # water fraction
296
297        nc_vegetation_type = self.nc_file.createVariable(
298            'vegetation_type', 'i1', ('y', 'x'), fill_value=-127)
299        nc_vegetation_type.long_name = "vegetation type classification"
300        nc_vegetation_type.units = "1"
301        nc_vegetation_type[:, :] = nc_vegetation_type._FillValue
302
303        nc_water_type = self.nc_file.createVariable(
304            'water_type', 'i1', ('y', 'x'), fill_value=-127)
305        nc_water_type.long_name = "water type classification"
306        nc_water_type.units = "1"
307        nc_water_type[:, :] = nc_water_type._FillValue
308
309        # Pixel-based surface and sub-surface parameters
310        # ----------------------------------------------
311        # nc_albedo_pars = self.nc_file.createVariable(
312        #     'albedo_pars', 'f4', ('nalbedo_pars', 'y', 'x'), fill_value=-9999.0)
313        # nc_albedo_pars.long_name = 'albedo parameters'
314        # nc_albedo_pars[:, :, :] = nc_albedo_pars._FillValue
315        #
316        # nc_building_pars = self.nc_file.createVariable(
317        #     'building_pars', 'f4', ('nbuilding_pars', 'y', 'x'), fill_value=-9999.0)
318        # nc_building_pars.long_name = 'building parameters'
319        # nc_building_pars[:, :, :] = nc_building_pars._FillValue
320        #
321        # nc_pavement_pars = self.nc_file.createVariable(
322        #     'pavement_pars', 'f4', ('npavement_pars', 'y', 'x'), fill_value=-9999.0)
323        # nc_pavement_pars.long_name = "pavement parameters"
324        # nc_pavement_pars[:, :, :] = nc_pavement_pars._FillValue
325        #
326        # nc_pavement_subsurface_pars = self.nc_file.createVariable(
327        #     'pavement_subsurface_pars', 'i1',
328        #     ('npavement_subsurface_pars','zsoil', 'y', 'x'), fill_value=-9999.0)
329        # nc_pavement_subsurface_pars.long_name = "pavement subsurface parameters"
330        # nc_pavement_subsurface_pars[:, :, :, :] = nc_pavement_subsurface_pars._FillValue
331        #
332        # nc_soil_pars = self.nc_file.createVariable(
333        #     'soil_pars', 'f', ('nsoil_pars','zsoil', 'y', 'x'), fill_value=-9999.0)
334        # nc_soil_pars.long_name = "soil parameters"
335        # nc_soil_pars.lod = np.int32(2) # if set to 1, adjust dimensions of soil_pars
336        # nc_soil_pars[:, :, :, :] = nc_soil_pars._FillValue
337        #
338        # nc_vegetation_pars = self.nc_file.createVariable(
339        #     'vegetation_pars', 'f', ('nvegetation_pars', 'y', 'x'), fill_value=-9999.0)
340        # nc_vegetation_pars.long_name = 'vegetation parameters'
341        # nc_vegetation_pars[:, :, :] = nc_vegetation_pars._FillValue
342        #
343        # nc_water_pars = self.nc_file.createVariable(
344        #     'water_pars', 'i1', ('nwater_pars', 'y', 'x'), fill_value=-9999.0)
345        # nc_water_pars.long_name = "water parameters"
346        # nc_water_pars[:, :, :] = nc_water_pars._FillValue
347
348        # Vegetation parameters
349        # ---------------------
350        nc_lad = self.nc_file.createVariable(
351            'lad', 'f4', ('zlad', 'y', 'x'), fill_value=-9999.0)
352        nc_lad.long_name = "leaf area density"
353        nc_lad.units = "m2 m-3"
354        nc_lad[:, :, :] = nc_lad._FillValue
355
356        # nc_bad = self.nc_file.createVariable(
357        #     'bad', 'f4', ('zlad', 'y', 'x'), fill_value=-9999.0)
358        # nc_bad.long_name = "basal area density"
359        # nc_bad.units = "m2 m-3"
360        # nc_bad[:, :, :] = nc_bad._FillValue
361        #
362        # nc_root_area_dens_r = self.nc_file.createVariable(
363        #     'root_area_dens_s', 'f4', ('zsoil', 'y', 'x'), fill_value=-9999.0)
364        # nc_root_area_dens_r.long_name = 'root area density of resolved vegetation'
365        # nc_root_area_dens_r.units = '1'
366        # nc_root_area_dens_r[:, :, :] = nc_root_area_dens_r._FillValue
367        #
368        # nc_root_area_dens_s = self.nc_file.createVariable(
369        #     'root_area_dens_s', 'f4', ('zsoil', 'y', 'x'), fill_value=-9999.0)
370        # nc_root_area_dens_s.long_name = 'root area density of parameterized vegetation'
371        # nc_root_area_dens_s.units = '1'
372        # nc_root_area_dens_s[:, :, :] = nc_root_area_dens_s._FillValue
373        #
374        # nc_tree_id = self.nc_file.createVariable(
375        #     'tree_id', 'i4', ('y', 'x'), fill_value=-9999.0)
376        # nc_tree_id.long_name = "tree id"
377        # nc_tree_id.units = "1"
378        # nc_tree_id[:, :, :] = nc_tree_id._FillValue
379
380        # Set topography
381        # --------------
382        nc_building_id[:5, :5] = 1
383        nc_building_id[-5:, :5] = 2
384        nc_building_id[:5, -5:] = 3
385        nc_building_id[-5:, -5:] = 4
386
387        nc_buildings_2d[:, :] = np.where(
388            nc_building_id[:, :] > nc_building_id._FillValue,
389            nc_building_id[:, :] * 10.0,
390            nc_buildings_2d[:, :])
391
392        for k in range(self.nz+1):
393            nc_buildings_3d[k, :, :] = np.where(nc_buildings_2d[:, :] > self.z[k], 1, 0)
394
395        nc_zt[:, :] = 4.0
396
397        # Set surface types
398        # -----------------
399        nc_building_type[:, :] = np.where(
400            nc_building_id[:, :] > nc_building_id._FillValue,
401            nc_building_id[:, :] * 1,
402            nc_building_type[:, :])
403
404        nc_pavement_type[9:11, :5] = 1
405        nc_pavement_type[:, 7:13] = 2
406
407        nc_vegetation_type[:, 5:7] = 3
408        nc_vegetation_type[:, 13:15] = 3
409        nc_vegetation_type[5:15, 15:] = 3
410        nc_vegetation_type[7:9, 15:17] = 1
411        nc_vegetation_type[11:13, 15:17] = 1
412
413        nc_water_type[5:9, :5] = 2
414        nc_water_type[11:15, :5] = 2
415
416        nc_soil_type[:, :] = np.where(
417            nc_vegetation_type[:, :] > nc_vegetation_type._FillValue,
418            2,
419            nc_soil_type[:, :])
420        nc_soil_type[:, :] = np.where(
421            nc_pavement_type[:, :] > nc_pavement_type._FillValue,
422            2,
423            nc_soil_type[:, :])
424
425        nc_street_type[:, :] = np.where(nc_pavement_type[:, :] == 1, 11, nc_street_type[:, :])
426        nc_street_type[:, :] = np.where(nc_pavement_type[:, :] == 2, 13, nc_street_type[:, :])
427
428        nc_surface_fraction[0, :, :] = np.where(
429            nc_building_id[:, :] > nc_building_id._FillValue,
430            nc_surface_fraction[0, :, :],
431            0)
432        nc_surface_fraction[2, :, :] = nc_surface_fraction[1, :, :] = nc_surface_fraction[0, :, :]
433        nc_surface_fraction[0, :, :] = np.where(
434            nc_vegetation_type[:, :] > nc_vegetation_type._FillValue,
435            1,
436            nc_surface_fraction[0, :, :])
437        nc_surface_fraction[1, :, :] = np.where(
438            nc_pavement_type[:, :] > nc_pavement_type._FillValue,
439            1,
440            nc_surface_fraction[1, :, :])
441        nc_surface_fraction[2, :, :] = np.where(
442            nc_water_type[:, :] > nc_water_type._FillValue,
443            1,
444            nc_surface_fraction[2, :, :])
445
446        # Set vegetation
447        # --------------
448        lad_profile = [0.0, 0.01070122, 0.1070122, 0.3130108, 0.3879193, 0.1712195]
449        for k in range(len(self.zlad)):
450            nc_lad[k, 7:9, 15:17] = lad_profile[k]
451            nc_lad[k, 11:13, 15:17] = lad_profile[k]
452
[4072]453    def finalize(self):
[4796]454        """Close file."""
[4072]455        print("Closing file...")
[4796]456
[4072]457        self.nc_file.close()
458
459
460if __name__ == '__main__':
461    driver = StaticDriver()
462    driver.write_global_attributes()
463    driver.define_dimensions()
[4796]464    driver.define_variables()
[4072]465    driver.finalize()
Note: See TracBrowser for help on using the repository browser.