Changeset 3629 for palm/trunk/SCRIPTS


Ignore:
Timestamp:
Dec 13, 2018 12:18:54 PM (5 years ago)
Author:
maronga
Message:

palm_csd improvements

Location:
palm/trunk/SCRIPTS
Files:
1 added
2 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SCRIPTS/palm_csd

    r3567 r3629  
    2525# -----------------
    2626# $Id$
     27# Added canopy generator calls. Some improvements
     28#
     29# 3567 2018-11-27 13:59:21Z maronga
    2730# Initial revisions
    28 #
    29 #
    30 #
    31 #
    3231#
    3332# Description:
     
    3837# @Author Bjoern Maronga (maronga@muk.uni-hannover.de)
    3938#
    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)
    4439# @todo Make input files optional
    4540# @todo Allow for ASCII input of terrain height and building height
     
    4944from palm_csd_files.palm_csd_netcdf_interface import *
    5045from palm_csd_files.palm_csd_tools import *
     46from palm_csd_files.palm_csd_canopy_generator import *
    5147import numpy as np
    5248
     
    7167      print ("Error. No configuration file " + input_config + " found.")
    7268      raise SystemExit   
    73    else:
    74       print(os.path.isfile(input_config))
    7569
    7670   config.read(input_config)
     
    7973#  Definition of settings
    8074   global settings_filename_out
    81    global settings_lai_season
    8275   global settings_bridge_width
     76   global settings_lai_roof_intensive
     77   global settings_lai_roof_extensive
     78   global settings_season
     79   global settings_lai_low_default
     80   global settings_lai_high_default
     81   global settings_patch_height_default
     82   global settings_lai_alpha
     83   global settings_lai_beta
    8384   global ndomains
    8485
     
    113114   global domain_dz
    114115   global domain_3d
    115    global domain_hv
    116    global domain_cg
     116   global domain_high_vegetation
    117117   global domain_ip
    118118   global domain_za
    119119   global domain_parent
     120   global domain_green_roofs
     121   global domain_street_trees
     122   global domain_canopy_patches
    120123
    121124#  Definition of input data parameters
     
    137140   global input_file_building_type
    138141   global input_file_building_type
     142   global input_file_lai 
    139143   global input_file_vegetation_type
    140144   global input_file_vegetation_height
     
    144148   global input_file_street_crossings
    145149   global input_file_soil_type
    146 
     150   global input_file_vegetation_on_roofs
     151   global input_file_tree_crown_diameter
     152   global input_file_tree_height
     153   global input_file_tree_trunk_diameter
     154   global input_file_tree_type
     155   global input_file_patch_height
    147156
    148157   global zt_all
     
    150159
    151160   settings_filename_out = "default"
    152    settings_lai_season = "summer"
    153161   settings_bridge_width = 3.0
     162   settings_lai_roof_intensive = 0.5
     163   settings_lai_roof_extensive = 1.0
     164   settings_season = "summer"
     165   settings_lai_high_default = 6.0
     166   settings_lai_low_default = 1.0
     167   settings_patch_height_default = 10.0
     168   settings_lai_alpha = 5.0
     169   settings_lai_beta = 3.0
    154170   ndomains = 0
    155171
     
    181197   domain_dz = []
    182198   domain_3d = []
    183    domain_hv = []
    184    domain_cg = []
     199   domain_high_vegetation = []
    185200   domain_ip = []
    186201   domain_za = []
    187202   domain_parent = []
     203   domain_green_roofs = []
     204   domain_street_trees = []
     205   domain_canopy_patches = []
    188206
    189207   zt_min = 0.0
     
    206224   input_file_bridges_id = []   
    207225   input_file_building_type = []
     226   input_file_lai = []
    208227   input_file_vegetation_type = []
    209228   input_file_vegetation_height = []   
     
    213232   input_file_street_crossings = []     
    214233   input_file_soil_type = []
    215 
     234   input_file_vegetation_on_roofs = []
     235   input_file_tree_crown_diameter = []
     236   input_file_tree_height = []
     237   input_file_tree_trunk_diameter = []
     238   input_file_tree_type = []
     239   input_file_patch_height = []
    216240   
    217241# Load all user parameters from config file
     
    240264      if ( read_tmp == 'settings' ):
    241265         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  
     266         settings_lai_roof_intensive = config.get(read_tmp, 'lai_roof_intensive')
     267         settings_lai_roof_extensive = config.get(read_tmp, 'lai_roof_extensive')
     268         settings_bridge_width = float(config.get(read_tmp, 'bridge_width'))   
     269         settings_season = config.get(read_tmp, 'season')   
     270         settings_lai_high_default = float(config.get(read_tmp, 'lai_high_vegetation_default'))
     271         settings_lai_low_default = float(config.get(read_tmp, 'lai_low_vegetation_default'))         
     272         settings_patch_height_default = float(config.get(read_tmp, 'patch_height_default')) 
     273         settings_lai_alpha = float(config.get(read_tmp, 'lai_alpha'))
     274         settings_lai_beta = float(config.get(read_tmp, 'lai_beta'))
     275
    245276      if ( read_tmp.split("_")[0] == 'domain' ):
    246277         ndomains = ndomains + 1
     
    251282         domain_dz.append(float(config.get(read_tmp, 'dz')))       
    252283         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'))       
     284         domain_high_vegetation.append(config.getboolean(read_tmp, 'allow_high_vegetation'))
     285         domain_canopy_patches.append(config.getboolean(read_tmp, 'generate_vegetation_patches'))       
    255286         domain_ip.append(config.getboolean(read_tmp, 'interpolate_terrain'))
    256287         domain_za.append(config.getboolean(read_tmp, 'use_palm_z_axis'))
     
    265296         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'))))
    266297         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 
     298         domain_green_roofs.append(config.getboolean(read_tmp, 'vegetation_on_roofs'))
     299         domain_street_trees.append(config.getboolean(read_tmp, 'street_trees'))
     300         
    268301      if ( read_tmp.split("_")[0] == 'input' ):
    269302         input_names.append(read_tmp.split("_")[1])
     
    280313         input_file_building_id.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_building_id')) 
    281314         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')) 
     315         input_file_building_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_building_type'))
     316         input_file_lai.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_lai'))         
    283317         input_file_vegetation_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_vegetation_type'))
    284318         input_file_vegetation_height.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_vegetation_height'))         
     
    286320         input_file_water_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_water_type'))
    287321         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')) 
     322         input_file_street_crossings.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_street_crossings'))
     323         input_file_patch_height.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_patch_height'))
     324         input_file_tree_crown_diameter.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_tree_crown_diameter')) 
     325         input_file_tree_height.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_tree_height')) 
     326         input_file_tree_trunk_diameter.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_tree_trunk_diameter')) 
     327         input_file_tree_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_tree_type')) 
     328         input_file_vegetation_on_roofs.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_vegetation_on_roofs'))           
    289329         #input_file_soil_type.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_soil_type'))   
    290330   return 0
     
    317357   "street_crossings": "b",     
    318358   "soil_type": "b",
    319    "surface_fraction": "f4"
     359   "surface_fraction": "f4",
     360   "building_pars": "f4",
     361   "vegetation_pars": "f4",
     362   "tree_data": "f4",
     363   "tree_type": "b",
     364   "nbuilding_pars": "i",
     365   "nvegetation_pars": "i",
     366   "zlad": "f4"
    320367   }
    321368
     
    340387   "street_crossings": np.byte(-127),   
    341388   "soil_type": np.byte(-127),
    342    "surface_fraction": float(-9999.0)
     389   "surface_fraction": float(-9999.0),
     390   "building_pars": float(-9999.0),
     391   "vegetation_pars": float(-9999.0),
     392   "tree_data": float(-9999.0),
     393   "tree_type": np.byte(-127)
    343394   }
    344395
     
    363414   "street_crossings": 0,   
    364415   "soil_type": 1,
    365    "surface_fraction": float(0.0)
     416   "surface_fraction": float(0.0),
     417   "buildings_pars": float(-9999.0),
     418   "tree_data": float(-9999.0),
     419   "tree_type": 0,
     420   "vegetation_pars": float(-9999.0)   
    366421   }
    367422 
     
    606661
    607662
    608       del bridges_2d, bridges_id, building_id
     663      del bridges_2d, bridges_id, building_id, buildings_2d
    609664
    610665
     
    684739   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])   
    685740         
    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)
     741   vegetation_type = np.where((vegetation_height != fillvalues["vegetation_height"]) & (vegetation_height == 0.0) & ((vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18)), 3, vegetation_type)
     742   vegetation_height = np.where((vegetation_height != fillvalues["vegetation_height"]) & (vegetation_height == 0.0) & ((vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18)), fillvalues["vegetation_height"],vegetation_height)
    688743
    689744
     
    720775   del soil_type
    721776
    722 
    723 
    724 
     777   del x
     778   del y
    725779
    726780#  pixels with bridges get building_type = 7 = bridge. This does not change the _type setting for the under-bridge
     
    737791         del bridges_2d
    738792
    739 # Read/Write street type and street crossings
     793# Read/write street type and street crossings
    740794for i in range(0,ndomains): 
    741795
     
    763817   nc_write_attribute(filename[i], 'street_crossings', 'grid_mapping', 'E_UTM N_UTM lon lat')   
    764818   del street_crossings
     819
     820
     821# Read/write vegetation on roofs
     822for i in range(0,ndomains): 
     823   if domain_green_roofs[i]:
     824      green_roofs = nc_read_from_file_2d(input_file_vegetation_on_roofs[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])   
     825      buildings_2d = nc_read_from_file_2d_all(filename[i], 'buildings_2d')
     826
     827
     828      x = nc_read_from_file_2d_all(filename[i], 'x')
     829      y = nc_read_from_file_2d_all(filename[i], 'y') 
     830      nbuilding_pars = np.arange(0,46)
     831      building_pars = np.ones((len(nbuilding_pars),len(y),len(x)))
     832      building_pars[:,:,:] = fillvalues["building_pars"]
     833
     834#     assign green fraction on roofs
     835      building_pars[3,:,:] = np.where( (buildings_2d != fillvalues["buildings_2d"] ) & ( green_roofs != fillvalues["building_pars"] ),1.0,fillvalues["building_pars"])
     836
     837#     assign leaf area index for vegetation on roofs
     838      building_pars[4,:,:] = np.where( ( buildings_2d != fillvalues["buildings_2d"] ) & ( green_roofs == 1.0 ),settings_lai_roof_intensive,fillvalues["building_pars"])
     839      building_pars[4,:,:] = np.where( ( buildings_2d != fillvalues["buildings_2d"] ) & ( green_roofs == 2.0 ),settings_lai_roof_extensive,building_pars[4,:,:])
     840   
     841   
     842      nc_write_dimension(filename[i], 'nbuilding_pars', nbuilding_pars, datatypes["nbuilding_pars"])
     843      nc_write_to_file_3d(filename[i], 'building_pars', building_pars, datatypes["building_pars"],'nbuilding_pars','y','x',fillvalues["building_pars"])
     844      nc_write_attribute(filename[i], 'building_pars', 'long_name', 'building_pars')
     845      nc_write_attribute(filename[i], 'building_pars', 'units', '')
     846      nc_write_attribute(filename[i], 'building_pars', 'res_orig', domain_px[i]) 
     847      nc_write_attribute(filename[i], 'building_pars', 'coordinates', 'E_UTM N_UTM lon lat')
     848      nc_write_attribute(filename[i], 'building_pars', 'grid_mapping', 'E_UTM N_UTM lon lat')   
     849         
     850      del building_pars, buildings_2d, x, y
     851
     852
     853# Read tree data and create LAD and BAD arrays using the canopy generator
     854for i in range(0,ndomains): 
     855   if domain_street_trees[i]:
     856      lai = nc_read_from_file_2d(input_file_lai[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])   
     857     
     858     
     859#     Save lai data as default for low and high vegetation
     860      lai_low = lai
     861      lai_high = lai
     862     
     863      x = nc_read_from_file_2d_all(filename[i], 'x')
     864      y = nc_read_from_file_2d_all(filename[i], 'y') 
     865      nvegetation_pars = np.arange(0,11)
     866      vegetation_pars = np.ones((len(nvegetation_pars),len(y),len(x)))
     867      vegetation_pars[:,:,:] = fillvalues["vegetation_pars"] 
     868     
     869      vegetation_pars[1,:,:] = lai
     870     
     871
     872#     Write out first version of LAI. Will later be overwritten.
     873      nc_write_dimension(filename[i], 'nvegetation_pars', nvegetation_pars, datatypes["nvegetation_pars"])
     874      nc_write_to_file_3d(filename[i], 'vegetation_pars', vegetation_pars, datatypes["vegetation_pars"],'nvegetation_pars','y','x',fillvalues["vegetation_pars"])
     875      nc_write_attribute(filename[i], 'vegetation_pars', 'long_name', 'vegetation_pars')
     876      nc_write_attribute(filename[i], 'vegetation_pars', 'units', '')
     877      nc_write_attribute(filename[i], 'vegetation_pars', 'res_orig', domain_px[i]) 
     878      nc_write_attribute(filename[i], 'vegetation_pars', 'coordinates', 'E_UTM N_UTM lon lat')
     879      nc_write_attribute(filename[i], 'vegetation_pars', 'grid_mapping', 'E_UTM N_UTM lon lat')   
     880
     881
     882#     Read all tree parameters from file
     883      tree_height = nc_read_from_file_2d(input_file_tree_height[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])   
     884      tree_crown_diameter = nc_read_from_file_2d(input_file_tree_crown_diameter[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])   
     885      tree_trunk_diameter = nc_read_from_file_2d(input_file_tree_trunk_diameter[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])         
     886      tree_type = nc_read_from_file_2d(input_file_tree_type[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])
     887      patch_height = nc_read_from_file_2d(input_file_patch_height[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])
     888     
     889#     Remove missing values from the data. Reasonable values will be set by the tree generator   
     890      tree_height = np.where( (tree_height == 0.0) | (tree_height == -1.0) ,fillvalues["tree_data"],tree_height)
     891      tree_crown_diameter = np.where( (tree_crown_diameter == 0.0) | (tree_crown_diameter == -1.0) ,fillvalues["tree_data"],tree_crown_diameter)
     892      tree_trunk_diameter = np.where( (tree_trunk_diameter == 0.0) | (tree_trunk_diameter == -1.0) ,fillvalues["tree_data"],tree_trunk_diameter)
     893      tree_type = np.where( (tree_type == 0.0) | (tree_type == -1.0) ,fillvalues["tree_data"],tree_type)
     894      patch_height = np.where( patch_height == -1.0 ,fillvalues["tree_data"],patch_height)     
     895
     896#     Convert trunk diameter from cm to m
     897      tree_trunk_diameter = np.where(tree_trunk_diameter != fillvalues["tree_data"], tree_trunk_diameter * 0.01,tree_trunk_diameter)
     898
     899
     900#     Temporarily change missing value for tree_type
     901      tree_type = np.where( (tree_type == fillvalues["tree_type"]),fillvalues["tree_data"],tree_type)   
     902
     903#     Compare patch height array with vegetation type and correct accordingly
     904      vegetation_type = nc_read_from_file_2d_all(filename[i], 'vegetation_type')
     905
     906
     907#     For zero-height patches, set vegetation_type to short grass and remove these pixels from the patch height array
     908      vegetation_type = np.where( (patch_height == 0.0) & ( (vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18) ),3,vegetation_type)
     909      patch_type = np.where( (patch_height == 0.0) & ( (vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18) ),fillvalues["tree_data"],patch_height)     
     910      nc_overwrite_to_file_2d(filename[i], 'vegetation_type', vegetation_type) 
     911
     912      max_tree_height = max(tree_height.flatten())
     913      max_patch_height = max(patch_height.flatten())
     914     
     915      if ( (max_tree_height != fillvalues["tree_data"]) | (max_patch_height == fillvalues["tree_data"]) ):
     916           
     917         lad, bad, tree_ids, zlad = generate_single_tree_lad(x,y,domain_dz[i],max_tree_height,max_patch_height,tree_type,tree_height,tree_crown_diameter,tree_trunk_diameter,lai,settings_season,fillvalues["tree_data"])
     918 
     919 
     920#        Remove LAD volumes that are inside buildings
     921         buildings_2d = nc_read_from_file_2d_all(filename[i], 'buildings_2d')
     922         for k in range(0,len(zlad)-1):
     923 
     924            lad[k,:,:] = np.where(buildings_2d == fillvalues["buildings_2d"],lad[k,:,:],fillvalues["tree_data"])
     925            bad[k,:,:] = np.where(buildings_2d == fillvalues["buildings_2d"],bad[k,:,:],fillvalues["tree_data"])
     926            tree_ids[k,:,:] = np.where(buildings_2d == fillvalues["buildings_2d"],tree_ids[k,:,:],fillvalues["tree_data"])     
     927 
     928 
     929         nc_write_dimension(filename[i], 'zlad', zlad, datatypes["tree_data"])
     930         nc_write_to_file_3d(filename[i], 'lad', lad, datatypes["tree_data"],'zlad','y','x',fillvalues["tree_data"])
     931         nc_write_attribute(filename[i], 'lad', 'long_name', 'leaf area density')
     932         nc_write_attribute(filename[i], 'lad', 'units', '')
     933         nc_write_attribute(filename[i], 'lad', 'res_orig', domain_px[i]) 
     934         nc_write_attribute(filename[i], 'lad', 'coordinates', 'E_UTM N_UTM lon lat')
     935         nc_write_attribute(filename[i], 'lad', 'grid_mapping', 'E_UTM N_UTM lon lat')
     936 
     937         nc_write_to_file_3d(filename[i], 'bad', bad, datatypes["tree_data"],'zlad','y','x',fillvalues["tree_data"])
     938         nc_write_attribute(filename[i], 'bad', 'long_name', 'basal area density')
     939         nc_write_attribute(filename[i], 'bad', 'units', '')
     940         nc_write_attribute(filename[i], 'bad', 'res_orig', domain_px[i]) 
     941         nc_write_attribute(filename[i], 'bad', 'coordinates', 'E_UTM N_UTM lon lat')
     942         nc_write_attribute(filename[i], 'bad', 'grid_mapping', 'E_UTM N_UTM lon lat')
     943           
     944         nc_write_to_file_3d(filename[i], 'tree_id', tree_ids, datatypes["tree_data"],'zlad','y','x',fillvalues["tree_data"])
     945         nc_write_attribute(filename[i], 'tree_id', 'long_name', 'tree id')
     946         nc_write_attribute(filename[i], 'tree_id', 'units', '')
     947         nc_write_attribute(filename[i], 'tree_id', 'res_orig', domain_px[i]) 
     948         nc_write_attribute(filename[i], 'tree_id', 'coordinates', 'E_UTM N_UTM lon lat')
     949         nc_write_attribute(filename[i], 'tree_id', 'grid_mapping', 'E_UTM N_UTM lon lat')           
     950           
     951         del buildings_2d, lai, lad, bad, tree_ids, zlad
     952
     953      del vegetation_pars, tree_height, tree_crown_diameter, tree_trunk_diameter, tree_type, patch_height, vegetation_type, x, y
     954
     955
     956# Create vegetation patches for locations with high vegetation type
     957for i in range(0,ndomains): 
     958   if domain_canopy_patches[i]:
     959
     960#     Load vegetation_type and lad array (at level z = 0) for re-processing
     961      vegetation_type = nc_read_from_file_2d_all(filename[i], 'vegetation_type')
     962      lad = nc_read_from_file_3d_all(filename[i], 'lad')
     963      patch_height = nc_read_from_file_2d(input_file_patch_height[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])
     964      vegetation_pars = nc_read_from_file_3d_all(filename[i], 'vegetation_pars')
     965      lai = vegetation_pars[1,:,:]
     966
     967#     For all pixels where the LAD array does not contains a street tree and for which a high vegetation was set and for which a patch height is either missing of larger than one grid spacing
     968#     the lai_high is reset to a low vegetation value
     969      lai_high = np.where( (lad[0,:,:] == fillvalues["tree_data"]) & ( ( (vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18) ) & ( (patch_height == fillvalues["tree_data"]) | (patch_height >= domain_dz[i])) ),settings_lai_low_default,fillvalues["tree_data"])
     970
     971#     For all pixels where lai_high is set but no lai is available, set the default high vegetation lai TODO: check if this makes sense bacauce lai_high = lai anway?!
     972      lai_high = np.where( (lai_high != fillvalues["tree_data"]) & (lai == fillvalues["tree_data"]), settings_lai_high_default, lai_high)
     973
     974#     Define a patch height whereever it is missing, but where a high vegetation LAI was set
     975      patch_height = np.where ( (lai_high != fillvalues["tree_data"]) & (patch_height == fillvalues["tree_data"]), settings_patch_height_default, patch_height)
     976
     977#     Remove pixels where street trees were already set
     978      patch_height = np.where ( (lad[0,:,:] != fillvalues["tree_data"]), fillvalues["tree_data"], patch_height)
     979
     980#     For missing LAI values, set either the high vegetation default or the low vegetation default   
     981      lai_high = np.where( (patch_height > 2.0) & (lai_high == fillvalues["tree_data"]),settings_lai_high_default,lai_high)
     982      lai_high = np.where( (patch_height <= 2.0) & (lai_high == fillvalues["tree_data"]),settings_lai_low_default,lai_high)
     983   
     984   
     985      if ( max(patch_height.flatten()) >= (2.0 * domain_dz[i]) ):
     986         print("    start calculating LAD (this might take some time)")
     987         lad_patch, patch_nz, status = process_patch(domain_dz[i],patch_height,lai_high,settings_lai_alpha,settings_lai_beta)
     988
     989         lad[0:patch_nz+1,:,:] = np.where( (lad[0:patch_nz+1,:,:] == fillvalues["tree_data"]),lad_patch[0:patch_nz+1,:,:], lad[0:patch_nz+1,:,:] )
     990
     991#     Remove high vegetation wherever it is replaced by a leaf area density. This should effectively remove all high vegetation pixels
     992      vegetation_type = np.where((lad[0,:,:] != fillvalues["tree_data"]) & (vegetation_type != fillvalues["vegetation_type"]),3,vegetation_type)
     993 
     994 
     995#     If desired, remove all high vegetation. TODO: check if this is still necessary
     996      if domain_high_vegetation[i]:
     997         vegetation_type = np.where((vegetation_type != fillvalues["vegetation_type"]) & ( (vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18) ),3,vegetation_type)   
     998
     999 
     1000#     Remove LAI from pixels where a LAD was set
     1001      lai_low = np.where( (lad[0,:,:] == fillvalues["tree_data"]), lai,fillvalues["tree_data"])
     1002   
     1003#     Fill low vegetation pixels without LAI set with default value
     1004      lai_low = np.where( (lai == fillvalues["tree_data"]) & (vegetation_type != fillvalues["vegetation_type"] ), settings_lai_low_default, lai)
     1005
     1006   
     1007#     Overwrite lai in vegetation_parameters 
     1008
     1009      vegetation_pars[1,:,:] = lai_low
     1010      nc_overwrite_to_file_3d(filename[i], 'vegetation_pars', vegetation_pars)
     1011
     1012#     Overwrite lad array
     1013      nc_overwrite_to_file_3d(filename[i], 'lad', lad)
     1014     
     1015#     Overwrite vegetation_type
     1016      nc_overwrite_to_file_2d(filename[i], 'vegetation_type', vegetation_type)   
     1017     
     1018     
     1019      del vegetation_type, lad, lai, patch_height, vegetation_pars
  • palm/trunk/SCRIPTS/palm_csd_files/palm_csd_netcdf_interface.py

    r3567 r3629  
    2525# -----------------
    2626# $Id$
     27# Some new routines
     28#
     29# 3567 2018-11-27 13:59:21Z maronga
    2730# Initial revisions
    28 #
    29 #
    30 #
    31 #
    3231#
    3332# Description:
     
    7675   nc_file = Dataset(filename, "r+", format="NETCDF4")
    7776   tmp_array = np.array(nc_file.variables[varname][:][:], dtype=type(nc_file.variables[varname]))
     77
     78   return tmp_array
     79
     80def nc_read_from_file_3d_all(filename, varname):
     81
     82
     83   import numpy as np
     84   import sys
     85   
     86   try:
     87      f = open(filename)
     88      f.close()
     89#      print("Load: " + filename + ".")
     90   except FileNotFoundError:
     91      print("Error: " + filename + ". No such file. Aborting...")
     92      sys.exit(1)
     93   
     94   nc_file = Dataset(filename, "r+", format="NETCDF4")
     95   tmp_array = np.array(nc_file.variables[varname][:][:][:], dtype=type(nc_file.variables[varname]))
    7896
    7997   return tmp_array
     
    212230   return 0
    213231
     232def nc_overwrite_to_file_3d(filename,varname,array):
     233
     234   try:
     235      f =  Dataset(filename, "a", format="NETCDF4")
     236      #print("Opened: " + filename + ".")
     237   except FileNotFoundError:
     238      print("Error. Could not open file: " + filename + ". Aborting...")
     239      sys.exit(1)
     240
     241   print("Writing array " + varname + " to file...")
     242   
     243   temp = f.variables[varname]
     244   temp[:,:,:] = array
     245   
     246   f.close()
     247
     248   return 0
     249
    214250def nc_write_to_file_3d(filename,varname,array,datatype,dimname1,dimname2,dimname3,fillvalue):
    215251
Note: See TracChangeset for help on using the changeset viewer.