Home


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

Changeset 3668 for palm/trunk/SCRIPTS


Ignore:
Timestamp:
Jan 14, 2019 12:49:24 PM (6 years ago)
Author:
maronga
Message:

removed most_methods circular and lookup. added improved version of palm_csd

Location:
palm/trunk/SCRIPTS
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SCRIPTS/.csd.config.default

    r3567 r3668  
    2525# -----------------
    2626# $Id$
     27# Updated configuration
     28#
     29# 3567 2018-11-27 13:59:21Z maronga
    2730# Initial revisions
    28 #
    29 #
    30 #
    31 #
    3231#
    3332# Description:
    … …  
    5251institution = Institute of Meteorology and Climatology, Leibniz University Hannover
    5352palm_version = 6.0
    54 version = 1
    5553rotation_angle = 0.0
    5654references = text
    … …  
    5856
    5957[settings]
    60 filename_out = /ldata2/MOSAIK/showcase_python
    61 lai_season = spring
    62 lai_high_vegetation_default = 5.0
    63 lai_low_vegetation_default = 1.0
     58lai_roof_intensive = 1.5
     59lai_roof_extensive = 3.0
     60lai_high_vegetation_default = 6.0
     61lai_low_vegetation_default = 3.0
    6462lai_alpha = 5.0
    6563lai_beta = 3.0
    66 tree_height_default = 10.0
     64patch_height_default = 10.0
    6765bridge_width = 3.0
    6866debug_mode = False
    69 
     67season = "summer"
    7068
    7169[input_01]
    … …  
    7977file_lon = Berlin_CoordinatesLatLon_x_15m_DLR.nc
    8078file_zt = Berlin_terrain_height_15m_DLR.nc
    81 file_buildings_2d = Berlin_building_height_15m_LOD2_DLR.nc
     79file_buildings_2d = Berlin_building_height_15m_DLR.nc
    8280file_building_id = Berlin_building_id_15m_DLR.nc
    8381file_building_type = Berlin_building_type_15m_DLR.nc
    8482file_bridges_2d = Berlin_bridges_height_15m_DLR.nc
    8583file_bridges_id = Berlin_bridges_id_15m_DLR.nc
     84file_lai =  Berlin_leaf_area_index_15m_DLR_WANG_summer.nc
    8685file_vegetation_type = Berlin_vegetation_type_15m_DLR.nc
    8786file_vegetation_height = Berlin_vegetation_patch_height_15m_DLR.nc
    … …  
    9190file_street_type =  Berlin_street_type_15m_DLR.nc
    9291file_street_crossings =  Berlin_street_crossings_15m_DLR.nc
     92file_tree_height =  Berlin_trees_height_clean_15m_DLR.nc
     93file_tree_crown_diameter =  Berlin_trees_crown_clean_15m_DLR.nc
     94file_tree_trunk_diameter =  Berlin_trees_trunk_clean_15m.nc
     95file_tree_type =  Berlin_trees_type_15m_DLR.nc
     96file_patch_height =  Berlin_vegetation_patch_height_15m_DLR.nc
     97file_vegetation_on_roofs =  Berlin_vegetation_on_roofs_15m_DLR.nc
    9398
    9499[input_02]
    … …  
    102107file_lon = Berlin_CoordinatesLatLon_x_10m_DLR.nc
    103108file_zt = Berlin_terrain_height_10m_DLR.nc
    104 file_buildings_2d = Berlin_building_height_10m_LOD2_DLR.nc
     109file_buildings_2d = Berlin_building_height_10m_DLR.nc
    105110file_building_id = Berlin_building_id_10m_DLR.nc
    106111file_building_type = Berlin_building_type_10m_DLR.nc
    107112file_bridges_2d = Berlin_bridges_height_10m_DLR.nc
    108113file_bridges_id = Berlin_bridges_id_10m_DLR.nc
     114file_lai =  Berlin_leaf_area_index_10m_DLR_WANG_summer.nc
    109115file_vegetation_type = Berlin_vegetation_type_10m_DLR.nc
    110116file_vegetation_height = Berlin_vegetation_patch_height_10m_DLR.nc
    … …  
    114120file_street_type =  Berlin_street_type_10m_DLR.nc
    115121file_street_crossings =  Berlin_street_crossings_10m_DLR.nc
     122file_tree_height =  Berlin_trees_height_clean_10m_DLR.nc
     123file_tree_crown_diameter =  Berlin_trees_crown_10m_DLR.nc
     124file_tree_trunk_diameter =  Berlin_trees_trunk_clean_10m.nc
     125file_tree_type =  Berlin_trees_type_10m_DLR.nc
     126file_patch_height =  Berlin_vegetation_patch_height_10m_DLR.nc
     127file_vegetation_on_roofs =  Berlin_vegetation_on_roofs_10m_DLR.nc
    116128
    117129[input_03]
     130path = /ldata2/MOSAIK/Berlin_static_driver_data
     131pixel_size = 1.0
    118132file_x = Berlin_CoordinatesUTM_x_1m_DLR.nc
    119133file_y = Berlin_CoordinatesUTM_y_1m_DLR.nc
    … …  
    122136file_lat = Berlin_CoordinatesLatLon_y_1m_DLR.nc
    123137file_lon = Berlin_CoordinatesLatLon_x_1m_DLR.nc
    124 path = /ldata2/MOSAIK/Berlin_static_driver_data
    125 pixel_size = 1.0
    126138file_zt = Berlin_terrain_height_1m_DLR.nc
    127 file_buildings_2d = Berlin_building_height_1m_LOD2_DLR.nc
     139file_buildings_2d = Berlin_building_height_1m_DLR.nc
    128140file_building_id = Berlin_building_id_1m_DLR.nc
    129141file_bridges_2d = Berlin_bridges_height_1m_DLR.nc
    130142file_bridges_id = Berlin_bridges_id_1m_DLR.nc
    131143file_building_type = Berlin_building_type_1m_DLR.nc
     144file_lai =  Berlin_leaf_area_index_1m_DLR_WANG_summer.nc
    132145file_vegetation_type = Berlin_vegetation_type_1m_DLR.nc
    133146file_vegetation_height = Berlin_vegetation_patch_height_1m_DLR.nc
    … …  
    137150file_street_type =  Berlin_street_type_1m_DLR.nc
    138151file_street_crossings =  Berlin_street_crossings_1m_DLR.nc
     152file_tree_height =  Berlin_trees_height_clean_1m_DLR.nc
     153file_tree_crown_diameter =  Berlin_trees_crown_clean_1m_DLR.nc
     154file_tree_trunk_diameter =  Berlin_trees_trunk_clean_1m.nc
     155file_tree_type =  Berlin_trees_type_1m_DLR.nc
     156file_patch_height =  Berlin_vegetation_patch_height_1m_DLR.nc
     157file_vegetation_on_roofs =  Berlin_vegetation_on_roofs_1m_DLR.nc
     158
     159[input_04]
     160path = /ldata2/MOSAIK/Berlin_static_driver_data
     161pixel_size = 2.0
     162file_x = Berlin_CoordinatesUTM_x_2m_DLR.nc
     163file_y = Berlin_CoordinatesUTM_y_2m_DLR.nc
     164file_x_UTM = Berlin_CoordinatesUTM_x_2m_DLR.nc
     165file_y_UTM = Berlin_CoordinatesUTM_y_2m_DLR.nc
     166file_lat = Berlin_CoordinatesLatLon_y_2m_DLR.nc
     167file_lon = Berlin_CoordinatesLatLon_x_2m_DLR.nc
     168file_zt = Berlin_terrain_height_2m_DLR.nc
     169file_buildings_2d = Berlin_building_height_2m_DLR.nc
     170file_building_id = Berlin_building_id_2m_DLR.nc
     171file_bridges_2d = Berlin_bridges_height_2m_DLR.nc
     172file_bridges_id = Berlin_bridges_id_2m_DLR.nc
     173file_building_type = Berlin_building_type_2m_DLR.nc
     174file_lai =  Berlin_leaf_area_index_2m_DLR_WANG_summer.nc
     175file_vegetation_type = Berlin_vegetation_type_2m_DLR.nc
     176file_vegetation_height = Berlin_vegetation_patch_height_2m_DLR.nc
     177file_pavement_type = Berlin_pavement_type_2m_DLR.nc
     178file_water_type = Berlin_water_type_2m_DLR.nc
     179file_soil_type = Berlin_soil_type_2m_DLR.nc
     180file_street_type =  Berlin_street_type_2m_DLR.nc
     181file_street_crossings =  Berlin_street_crossings_2m_DLR.nc
     182file_tree_height =  Berlin_trees_height_clean_2m_DLR.nc
     183file_tree_crown_diameter = Berlin_trees_crown_clean_2m_DLR.nc
     184file_tree_trunk_diameter =  Berlin_trees_trunk_clean_2m.nc
     185file_tree_type =  Berlin_trees_type_2m_DLR.nc
     186file_patch_height =  Berlin_vegetation_patch_height_2m_DLR.nc
     187file_vegetation_on_roofs =  Berlin_vegetation_on_roofs_2m_DLR.nc
     188
    139189
    140190[output]
    141191path = /ldata2/MOSAIK/
    142 file_out = winter_iop1_test
     192file_out = showcase_berlin
    143193version = 1
    144194
    145195[domain_root]
    146196pixel_size = 15.0
    147 origin_x = 0.0
    148 origin_y = 0.0
    149 nx = 311
    150 ny = 257
     197origin_x = 0
     198origin_y = 0
     199nx = 3119
     200ny = 2573
    151201buildings_3d = False
    152202dz = 15.0
    153 allow_high_vegetation = True
     203allow_high_vegetation = False
    154204generate_vegetation_patches = True
    155 use_palm_z_axis= True
     205use_palm_z_axis= False
    156206interpolate_terrain = False
    157207domain_parent
     208vegetation_on_roofs = False
     209street_trees = True
    158210
    159211[domain_N02]
    … …  
    170222interpolate_terrain = True
    171223domain_parent = root
     224vegetation_on_roofs = False
     225street_trees = True

  • palm/trunk/SCRIPTS/palm_csd

    r3629 r3668  
    2525# -----------------
    2626# $Id$
     27# Various improvements and bugfixes
     28#
     29# 3629 2018-12-13 12:18:54Z maronga
    2730# Added canopy generator calls. Some improvements
    2831#
    … …  
    7275
    7376#  Definition of settings
    74    global settings_filename_out
     77
    7578   global settings_bridge_width
    7679   global settings_lai_roof_intensive
    … …  
    100103   global global_site
    101104   global global_source
    102    global global_version
     105
     106   global  path_out
     107   global filename_out
     108   global version_out
    103109
    104110
    … …  
    158164   global zt_min
    159165
    160    settings_filename_out = "default"
    161166   settings_bridge_width = 3.0
    162167   settings_lai_roof_intensive = 0.5
    … …  
    185190   global_site = ""
    186191   global_source = ""
    187    global_version = 1
    188 
     192
     193   path_out = ""
     194   version_out = 1
     195   filename_out = "default"
     196   
    189197   domain_names = []
    190198   domain_px = []
    … …  
    260268         global_site = config.get(read_tmp, 'site')
    261269         global_source = config.get(read_tmp, 'source')
    262          global_version = float(config.get(read_tmp, 'version'))
     270
    263271
    264272      if ( read_tmp == 'settings' ):
    265          settings_filename_out = config.get(read_tmp, 'filename_out')
    266273         settings_lai_roof_intensive = config.get(read_tmp, 'lai_roof_intensive')
    267274         settings_lai_roof_extensive = config.get(read_tmp, 'lai_roof_extensive')
    … …  
    273280         settings_lai_alpha = float(config.get(read_tmp, 'lai_alpha'))
    274281         settings_lai_beta = float(config.get(read_tmp, 'lai_beta'))
     282
     283      if ( read_tmp == 'output' ):
     284         path_out = config.get(read_tmp, 'path')         
     285         filename_out = config.get(read_tmp, 'file_out')
     286         version_out = float(config.get(read_tmp, 'version'))
    275287
    276288      if ( read_tmp.split("_")[0] == 'domain' ):
    … …  
    322334         input_file_street_crossings.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_street_crossings'))
    323335         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')) 
     336         
     337         tmp = config.get(read_tmp, 'file_tree_crown_diameter')
     338         if tmp is not None:
     339            input_file_tree_crown_diameter.append(config.get(read_tmp, 'path') + "/" + tmp) 
     340         else:
     341            input_file_tree_crown_diameter.append(None)
    325342         input_file_tree_height.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_tree_height')) 
    326343         input_file_tree_trunk_diameter.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_tree_trunk_diameter')) 
    … …  
    433450#  Calculate indices and input files
    434451   ii.append(input_px.index(domain_px[i]))
    435    filename.append(settings_filename_out + "_" + domain_names[i])
     452   filename.append(path_out + "/" + filename_out + "_" + domain_names[i])
    436453   if domain_parent[i] is not None:
    437454      ii_parent.append(domain_names.index(domain_parent[i]))
    … …  
    447464#  Create NetCDF output file and set global attributes
    448465   nc_create_file(filename[i])
    449    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)
     466   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,version_out)
    450467
    451468   
    … …  
    468485del zt_all[:]
    469486
     487print( "Shift terrain height by -" + str(zt_min))
    470488for i in range(0,ndomains):
    471489
    … …  
    475493   y = nc_read_from_file_1d(input_file_y[ii[i]], "y", domain_y0[i], domain_y1[i])
    476494
    477    print( "Shift terrain height by -" + str(zt_min))
    478495   zt = zt - zt_min
     496   
     497   nc_write_global_attribute(filename[i], 'origin_z', float(zt_min))
    479498
    480499#  If necessary, interpolate parent domain terrain height on child domain grid and blend the two
    … …  
    494513
    495514#     Interpolate array and bring to PALM grid of child domain
     515
    496516      zt_ip = interpolate_2d(zt_parent,tmp_x,tmp_y,x,y)
    497517      zt_ip = bring_to_palm_grid(zt_ip,x,y,domain_dz[ii_parent[i]])
    … …  
    580600 
    581601   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]) 
    582    building_type[building_type == 255] = fillvalues["building_type"]
     602   building_type[building_type >= 254] = fillvalues["building_type"]
    583603   building_type = np.where(building_type < 1,defaultvalues["building_type"],building_type)
    584604 
    … …  
    594614      building_type = np.where((building_type == fillvalues["building_type"]) & (buildings_2d != fillvalues["buildings_2d"]),defaultvalues["building_type"],building_type)
    595615      print("Data check #2 " + str(check_arrays_2(buildings_2d,building_type,fillvalues["buildings_2d"],fillvalues["building_type"])))
     616 
    596617 
    597618   nc_write_to_file_2d(filename[i], 'buildings_2d', buildings_2d, datatypes["buildings_2d"],'y','x',fillvalues["buildings_2d"])
    … …  
    606627   nc_write_attribute(filename[i], 'building_id', 'long_name', 'building id')
    607628   nc_write_attribute(filename[i], 'building_id', 'units', '')
    608    nc_write_attribute(filename[i], 'building_id', 'res_orig', domain_px[i])
     629   nc_write_attribute(filename[i], 'building_id', 'res _orig', domain_px[i])
    609630   nc_write_attribute(filename[i], 'building_id', 'coordinates', 'E_UTM N_UTM lon lat')
    610631   nc_write_attribute(filename[i], 'building_id', 'grid_mapping', 'E_UTM N_UTM lon lat')   
    … …  
    616637   nc_write_attribute(filename[i], 'building_type', 'coordinates', 'E_UTM N_UTM lon lat')
    617638   nc_write_attribute(filename[i], 'building_type', 'grid_mapping', 'E_UTM N_UTM lon lat')   
     639 
     640   print("out_1: " + str(np.any(building_type == -2)))
    618641 
    619642del buildings_2d
    … …  
    705728#  #6 Remove water for pixels with buildings
    706729   water_type = np.where((water_type != fillvalues["water_type"]) & (building_type != fillvalues["building_type"]),fillvalues["water_type"],water_type)
    707    
    708 
    709 #  #7 to be removed: set default soil type everywhere
    710    soil_type = np.where((vegetation_type != fillvalues["vegetation_type"]) | (pavement_type != fillvalues["pavement_type"]),defaultvalues["soil_type"],fillvalues["soil_type"])
    711  
     730 
     731   
     732#  Correct vegetation_type when a vegetation height is available and is indicative of low vegeetation
     733   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])   
     734         
     735   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)
     736   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)
    712737 
    713738#  Check for consistency and fill empty fields with default vegetation type
    … …  
    717742      vegetation_type = np.where(consistency_array == 0,defaultvalues["vegetation_type"],vegetation_type)
    718743      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"])
    719    
     744
     745#  #7 to be removed: set default soil type everywhere
     746   soil_type = np.where((vegetation_type != fillvalues["vegetation_type"]) | (pavement_type != fillvalues["pavement_type"]),defaultvalues["soil_type"],fillvalues["soil_type"])
     747       
     748   
     749#  Check for consistency and fill empty fields with default vegetation type
     750   consistency_array, test = check_consistency_3(vegetation_type,pavement_type,soil_type,fillvalues["vegetation_type"],fillvalues["pavement_type"],fillvalues["soil_type"])
     751       
    720752#  Create surface_fraction array
    721753   x = nc_read_from_file_2d_all(filename[i], 'x')
    … …  
    734766   nc_write_attribute(filename[i], 'surface_fraction', 'res_orig', domain_px[i])     
    735767   del surface_fraction   
    736    
    737  
    738 #  Correct vegetation_type when a vegetation height is available and is indicative of low vegeetation
    739    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])   
    740          
    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)
    743 
    744768
    745769   nc_write_to_file_2d(filename[i], 'vegetation_type', vegetation_type, datatypes["vegetation_type"],'y','x',fillvalues["vegetation_type"])     
    … …  
    779803
    780804#  pixels with bridges get building_type = 7 = bridge. This does not change the _type setting for the under-bridge
    781 #  area
     805#  area NOTE: when bridges are present the consistency check will fail at the moment
    782806   if domain_3d[i]:
    783807      if np.any(building_type != fillvalues["building_type"]):
    … …  
    853877# Read tree data and create LAD and BAD arrays using the canopy generator
    854878for i in range(0,ndomains): 
     879   lai = nc_read_from_file_2d(input_file_lai[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])   
     880     
     881   vegetation_type = nc_read_from_file_2d_all(filename[i], 'vegetation_type')
     882     
     883   lai = np.where(vegetation_type == fillvalues["vegetation_type"],fillvalues["vegetation_pars"],lai)
     884
     885     
     886   x = nc_read_from_file_2d_all(filename[i], 'x')
     887   y = nc_read_from_file_2d_all(filename[i], 'y') 
     888   nvegetation_pars = np.arange(0,12)
     889   vegetation_pars = np.ones((len(nvegetation_pars),len(y),len(x)))
     890   vegetation_pars[:,:,:] = fillvalues["vegetation_pars"] 
     891     
     892   vegetation_pars[1,:,:] = lai
     893     
     894
     895#  Write out first version of LAI. Will later be overwritten.
     896   nc_write_dimension(filename[i], 'nvegetation_pars', nvegetation_pars, datatypes["nvegetation_pars"])
     897   nc_write_to_file_3d(filename[i], 'vegetation_pars', vegetation_pars, datatypes["vegetation_pars"],'nvegetation_pars','y','x',fillvalues["vegetation_pars"])
     898   nc_write_attribute(filename[i], 'vegetation_pars', 'long_name', 'vegetation_pars')
     899   nc_write_attribute(filename[i], 'vegetation_pars', 'units', '')
     900   nc_write_attribute(filename[i], 'vegetation_pars', 'res_orig', domain_px[i]) 
     901   nc_write_attribute(filename[i], 'vegetation_pars', 'coordinates', 'E_UTM N_UTM lon lat')
     902   nc_write_attribute(filename[i], 'vegetation_pars', 'grid_mapping', 'E_UTM N_UTM lon lat')   
     903
     904   del lai, vegetation_pars, vegetation_type
     905
     906# Read tree data and create LAD and BAD arrays using the canopy generator
     907for i in range(0,ndomains): 
    855908   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      
     909
     910      vegetation_pars =  nc_read_from_file_2d_all(filename[i], 'vegetation_pars')
     911
     912      lai = np.copy(vegetation_pars[1,:,:])
     913
     914      x = nc_read_from_file_2d_all(filename[i], 'x')
     915      y = nc_read_from_file_2d_all(filename[i], 'y') 
     916
    859917#     Save lai data as default for low and high vegetation
    860918      lai_low = lai
    861919      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 
     920   
    882921#     Read all tree parameters from file
    883922      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])   
     923
     924      if (input_file_tree_crown_diameter[ii[i]] is not None):
     925         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])   
     926         tree_crown_diameter = np.where( (tree_crown_diameter == 0.0) | (tree_crown_diameter == -1.0) ,fillvalues["tree_data"],tree_crown_diameter)
     927      else:
     928         tree_crown_diameter = np.ones((len(y),len(x)))
     929         tree_crown_diameter[:,:] = fillvalues["tree_data"]
     930
     931
    885932      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])         
    886933      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])
    … …  
    889936#     Remove missing values from the data. Reasonable values will be set by the tree generator   
    890937      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)
    892938      tree_trunk_diameter = np.where( (tree_trunk_diameter == 0.0) | (tree_trunk_diameter == -1.0) ,fillvalues["tree_data"],tree_trunk_diameter)
    893939      tree_type = np.where( (tree_type == 0.0) | (tree_type == -1.0) ,fillvalues["tree_data"],tree_type)
    … …  
    907953#     For zero-height patches, set vegetation_type to short grass and remove these pixels from the patch height array
    908954      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) 
     955      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)   
     956     
    911957
    912958      max_tree_height = max(tree_height.flatten())
    … …  
    914960     
    915961      if ( (max_tree_height != fillvalues["tree_data"]) | (max_patch_height == fillvalues["tree_data"]) ):
    916            
     962
    917963         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"])
    918964 
    … …  
    925971            bad[k,:,:] = np.where(buildings_2d == fillvalues["buildings_2d"],bad[k,:,:],fillvalues["tree_data"])
    926972            tree_ids[k,:,:] = np.where(buildings_2d == fillvalues["buildings_2d"],tree_ids[k,:,:],fillvalues["tree_data"])     
    927  
     973
     974         del buildings_2d
    928975 
    929976         nc_write_dimension(filename[i], 'zlad', zlad, datatypes["tree_data"])
    … …  
    949996         nc_write_attribute(filename[i], 'tree_id', 'grid_mapping', 'E_UTM N_UTM lon lat')           
    950997           
    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
     998         del lai, lad, bad, tree_ids, zlad
     999
     1000      del vegetation_pars, tree_height, tree_crown_diameter, tree_trunk_diameter, tree_type, patch_height, x, y
    9541001
    9551002
    … …  
    9611008      vegetation_type = nc_read_from_file_2d_all(filename[i], 'vegetation_type')
    9621009      lad = nc_read_from_file_3d_all(filename[i], 'lad')
     1010      zlad = nc_read_from_file_1d_all(filename[i], 'zlad')
    9631011      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])
    9641012      vegetation_pars = nc_read_from_file_3d_all(filename[i], 'vegetation_pars')
    9651013      lai = vegetation_pars[1,:,:]
    9661014
    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
     1015
     1016#     Determine all pixels that do not already have and LAD but which are high vegetation to a dummy value of 1.0 and remove all other pixels
     1017      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])) ),1.0,fillvalues["tree_data"])
     1018
     1019#     Now, assign either the default LAI for high vegetation or the value stored in the LAI array.
     1020      lai_high = np.where( (lai_high != fillvalues["tree_data"]) & (lai == fillvalues["tree_data"]), settings_lai_high_default, lai)
     1021
     1022#     Define a patch height wherever it is missing, but where a high vegetation LAI was set
    9751023      patch_height = np.where ( (lai_high != fillvalues["tree_data"]) & (patch_height == fillvalues["tree_data"]), settings_patch_height_default, patch_height)
    9761024
    … …  
    9781026      patch_height = np.where ( (lad[0,:,:] != fillvalues["tree_data"]), fillvalues["tree_data"], patch_height)
    9791027
     1028#     Remove patch heights that have no lai_high value
     1029      patch_height = np.where ( (lai_high == fillvalues["tree_data"]), fillvalues["tree_data"], patch_height)
     1030
    9801031#     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    
     1032      lai_high = np.where( (patch_height > 2.0) & (patch_height != fillvalues["tree_data"]) & (lai_high == fillvalues["tree_data"]),settings_lai_high_default,lai_high)
     1033      lai_high = np.where( (patch_height <= 2.0) & (patch_height != fillvalues["tree_data"]) & (lai_high == fillvalues["tree_data"]),settings_lai_low_default,lai_high)
     1034         
    9851035      if ( max(patch_height.flatten()) >= (2.0 * domain_dz[i]) ):
    9861036         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)
     1037         lad_patch, patch_nz, status = process_patch(domain_dz[i],patch_height,max(zlad),lai_high,settings_lai_alpha,settings_lai_beta)
    9881038
    9891039         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,:,:] )
    … …  
    9911041#     Remove high vegetation wherever it is replaced by a leaf area density. This should effectively remove all high vegetation pixels
    9921042      vegetation_type = np.where((lad[0,:,:] != fillvalues["tree_data"]) & (vegetation_type != fillvalues["vegetation_type"]),3,vegetation_type)
    993  
    994  
     1043             
    9951044#     If desired, remove all high vegetation. TODO: check if this is still necessary
    9961045      if domain_high_vegetation[i]:
    … …  
    9981047
    9991048 
    1000 #     Remove LAI from pixels where a LAD was set
    1001       lai_low = np.where( (lad[0,:,:] == fillvalues["tree_data"]), lai,fillvalues["tree_data"])
     1049#     Set default low LAI for pixels with an LAD (short grass below trees)
     1050      lai_low = np.where( (lad[0,:,:] == fillvalues["tree_data"]), lai, settings_lai_low_default)
    10021051   
    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    
     1052#     Fill low vegetation pixels without LAI set or with LAI = 0 with default value
     1053      lai_low = np.where( ( (lai_low == fillvalues["tree_data"]) | (lai_low == 0.0) ) & (vegetation_type != fillvalues["vegetation_type"] ), settings_lai_low_default, lai_low)
     1054
     1055#     Remove lai for pixels that have no vegetation_type
     1056      lai_low = np.where( vegetation_type != fillvalues["vegetation_type"], lai_low, fillvalues["tree_data"])
     1057       
    10071058#     Overwrite lai in vegetation_parameters 
    1008 
    1009       vegetation_pars[1,:,:] = lai_low
     1059      vegetation_pars[1,:,:] = np.copy(lai_low)
    10101060      nc_overwrite_to_file_3d(filename[i], 'vegetation_pars', vegetation_pars)
    10111061
    10121062#     Overwrite lad array
    10131063      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
     1064           
     1065      nc_overwrite_to_file_2d(filename[i], 'vegetation_type', vegetation_type) 
     1066
     1067     
     1068      del vegetation_type, lad, lai, patch_height, vegetation_pars, zlad
     1069
     1070# Final consistency check
     1071for i in range(0,ndomains):
     1072   vegetation_type = nc_read_from_file_2d_all(filename[i], 'vegetation_type') 
     1073   pavement_type = nc_read_from_file_2d_all(filename[i], 'pavement_type') 
     1074   building_type = nc_read_from_file_2d_all(filename[i], 'building_type') 
     1075   water_type = nc_read_from_file_2d_all(filename[i], 'water_type')     
     1076   soil_type = nc_read_from_file_2d_all(filename[i], 'soil_type') 
     1077
     1078#  Check for consistency and fill empty fields with default vegetation type
     1079   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"])
     1080
     1081#  Check for consistency and fill empty fields with default vegetation type
     1082   consistency_array, test = check_consistency_3(vegetation_type,pavement_type,soil_type,fillvalues["vegetation_type"],fillvalues["pavement_type"],fillvalues["soil_type"])
     1083
     1084   surface_fraction = nc_read_from_file_3d_all(filename[i], 'surface_fraction')       
     1085   surface_fraction[0,:,:] = np.where(vegetation_type != fillvalues["vegetation_type"], 1.0, 0.0)
     1086   surface_fraction[1,:,:] = np.where(pavement_type != fillvalues["pavement_type"], 1.0, 0.0)
     1087   surface_fraction[2,:,:] = np.where(water_type != fillvalues["water_type"], 1.0, 0.0)
     1088   nc_overwrite_to_file_3d(filename[i], 'surface_fraction', surface_fraction)
     1089       
     1090   del vegetation_type, pavement_type, building_type, water_type, soil_type

  • palm/trunk/SCRIPTS/palm_csd_files/palm_csd_canopy_generator.py

    r3629 r3668  
    2525# -----------------
    2626# $Id$
     27# Various improvements and bugfixes
     28#
     29# 3629 2018-12-13 12:18:54Z maronga
    2730# Initial revision
    2831#
    … …  
    7275            if valid_pixels[j,i]:
    7376               tree_id_counter = tree_id_counter + 1
    74                print("   Processing tree No " +  str(tree_id_counter) + " ...", end="")
    75 
     77#               print("   Processing tree No " +  str(tree_id_counter) + " ...", end="")
    7678               lad_loc, bad_loc, x_loc, y_loc, z_loc, status = process_single_tree(dx,dz,tree_type[j,i],tree_height[j,i],tree_lai[j,i],tree_dia[j,i],trunk_dia[j,i],season)
    7779               if ( np.any(lad_loc) != fill ):
    … …  
    101103
    102104
    103                   if ( status == 0 ):
    104                      status_char = " ok."
    105                   else:
    106                      status_char = " skipped."
    107                   print(status_char)
     105#                  if ( status == 0 ):
     106#                     status_char = " ok."
     107#                  else:
     108#                     status_char = " skipped."
     109#                  print(status_char)
    108110
    109111               del lad_loc, x_loc, y_loc, z_loc, status
    … …  
    131133   default_trees = []
    132134                                                #1   #2   #3   #4    #5   #6   #7   #8     #9
    133    default_trees.append(tree("Default",         1.0, 1.0, 4.0, 12.0, 2.0, 0.1, 0.6, 0.025, 0.35))
    134    default_trees.append(tree("Abies",           3.0, 1.0, 4.0, 12.0, 2.0, 0.1, 0.6, 0.025, 0.80))
    135    default_trees.append(tree("Acer",            1.0, 1.0, 7.0, 12.0, 2.0, 0.1, 0.6, 0.025, 0.80))
    136    default_trees.append(tree("Aesculus",        1.0, 1.0, 7.0, 12.0, 2.0, 0.1, 0.6, 0.025, 1.00))
    137    default_trees.append(tree("Ailanthus",       1.0, 1.0, 8.5, 13.5, 2.0, 0.1, 0.6, 0.025, 1.30))
    138    default_trees.append(tree("Alnus",           3.0, 1.0, 6.0, 16.0, 2.0, 0.1, 0.6, 0.025, 1.20))
    139    default_trees.append(tree("Amelanchier",     1.0, 1.0, 3.0,  4.0, 2.0, 0.1, 0.6, 0.025, 1.20))
    140    default_trees.append(tree("Betula",          1.0, 1.0, 6.0, 14.0, 2.0, 0.1, 0.6, 0.025, 0.30))
    141    default_trees.append(tree("Buxus",           1.0, 1.0, 4.0,  4.0, 2.0, 0.1, 0.6, 0.025, 0.90))
    142    default_trees.append(tree("Calocedrus",      3.0, 1.0, 5.0, 10.0, 2.0, 0.1, 0.6, 0.025, 0.50))
    143    default_trees.append(tree("Caragana",        1.0, 1.0, 3.5,  6.0, 2.0, 0.1, 0.6, 0.025, 0.90))
    144    default_trees.append(tree("Carpinus",        1.0, 1.0, 6.0, 10.0, 2.0, 0.1, 0.6, 0.025, 0.70))
    145    default_trees.append(tree("Carya",           1.0, 1.0, 5.0, 17.0, 2.0, 0.1, 0.6, 0.025, 0.80))
    146    default_trees.append(tree("Castanea",        1.0, 1.0, 4.5,  7.0, 2.0, 0.1, 0.6, 0.025, 0.80))
    147    default_trees.append(tree("Catalpa",         1.0, 1.0, 5.5,  6.5, 2.0, 0.1, 0.6, 0.025, 0.70))
    148    default_trees.append(tree("Cedrus",          1.0, 1.0, 8.0, 13.0, 2.0, 0.1, 0.6, 0.025, 0.80))
    149    default_trees.append(tree("Celtis",          1.0, 1.0, 6.0,  9.0, 2.0, 0.1, 0.6, 0.025, 0.80))
    150    default_trees.append(tree("Cercidiphyllum",  1.0, 1.0, 3.0,  6.5, 2.0, 0.1, 0.6, 0.025, 0.80))
    151    default_trees.append(tree("Cercis",          1.0, 1.0, 2.5,  7.5, 2.0, 0.1, 0.6, 0.025, 0.90))
    152    default_trees.append(tree("Chamaecyparis",   5.0, 1.0, 3.5,  9.0, 2.0, 0.1, 0.6, 0.025, 0.70))
    153    default_trees.append(tree("Cladrastis",      1.0, 1.0, 5.0, 10.0, 2.0, 0.1, 0.6, 0.025, 0.80))
    154    default_trees.append(tree("Cornus",          1.0, 1.0, 4.5,  6.5, 2.0, 0.1, 0.6, 0.025, 1.20))
    155    default_trees.append(tree("Corylus",         1.0, 1.0, 5.0,  9.0, 2.0, 0.1, 0.6, 0.025, 0.40))
    156    default_trees.append(tree("Cotinus",         1.0, 1.0, 4.0,  4.0, 2.0, 0.1, 0.6, 0.025, 0.70))
    157    default_trees.append(tree("Crataegus",       3.0, 1.0, 3.5,  6.0, 2.0, 0.1, 0.6, 0.025, 1.40))
    158    default_trees.append(tree("Cryptomeria",     3.0, 1.0, 5.0, 10.0, 2.0, 0.1, 0.6, 0.025, 0.50))
    159    default_trees.append(tree("Cupressocyparis", 3.0, 1.0, 3.0,  8.0, 2.0, 0.1, 0.6, 0.025, 0.40))
    160    default_trees.append(tree("Cupressus",       3.0, 1.0, 5.0,  7.0, 2.0, 0.1, 0.6, 0.025, 0.40))
    161    default_trees.append(tree("Cydonia",         1.0, 1.0, 2.0,  3.0, 2.0, 0.1, 0.6, 0.025, 0.90))
    162    default_trees.append(tree("Davidia",         1.0, 1.0,10.0, 14.0, 2.0, 0.1, 0.6, 0.025, 0.40))
    163    default_trees.append(tree("Elaeagnus",       1.0, 1.0, 6.5,  6.0, 2.0, 0.1, 0.6, 0.025, 1.20))
    164    default_trees.append(tree("Euodia",          1.0, 1.0, 4.5,  6.0, 2.0, 0.1, 0.6, 0.025, 0.90))
    165    default_trees.append(tree("Euonymus",        1.0, 1.0, 4.5,  6.0, 2.0, 0.1, 0.6, 0.025, 0.60))
    166    default_trees.append(tree("Fagus",           1.0, 1.0,10.0, 12.5, 2.0, 0.1, 0.6, 0.025, 0.50))
    167    default_trees.append(tree("Fraxinus",        1.0, 1.0, 5.5, 10.5, 2.0, 0.1, 0.6, 0.025, 1.60))
    168    default_trees.append(tree("Ginkgo",          3.0, 1.0, 4.0,  8.5, 2.0, 0.1, 0.6, 0.025, 0.80))
    169    default_trees.append(tree("Gleditsia",       1.0, 1.0, 6.5, 10.5, 2.0, 0.1, 0.6, 0.025, 0.60))
    170    default_trees.append(tree("Gymnocladus",     1.0, 1.0, 5.5, 10.0, 2.0, 0.1, 0.6, 0.025, 0.80))
    171    default_trees.append(tree("Hippophae",       1.0, 1.0, 9.5,  8.5, 2.0, 0.1, 0.6, 0.025, 0.80))
    172    default_trees.append(tree("Ilex",            1.0, 1.0, 4.0,  7.5, 2.0, 0.1, 0.6, 0.025, 0.80))
    173    default_trees.append(tree("Juglans",         1.0, 1.0, 7.0,  9.0, 2.0, 0.1, 0.6, 0.025, 0.50))
    174    default_trees.append(tree("Juniperus",       5.0, 1.0, 3.0,  7.0, 2.0, 0.1, 0.6, 0.025, 0.90))
    175    default_trees.append(tree("Koelreuteria",    1.0, 1.0, 3.5,  5.5, 2.0, 0.1, 0.6, 0.025, 0.50))
    176    default_trees.append(tree("Laburnum",        1.0, 1.0, 3.0,  6.0, 2.0, 0.1, 0.6, 0.025, 0.60))
    177    default_trees.append(tree("Larix",           3.0, 1.0, 7.0, 16.5, 2.0, 0.1, 0.6, 0.025, 0.60))
    178    default_trees.append(tree("Ligustrum",       1.0, 1.0, 3.0,  6.0, 2.0, 0.1, 0.6, 0.025, 1.10))
    179    default_trees.append(tree("Liquidambar",     3.0, 1.0, 3.0,  7.0, 2.0, 0.1, 0.6, 0.025, 0.30))
    180    default_trees.append(tree("Liriodendron",    3.0, 1.0, 4.5,  9.5, 2.0, 0.1, 0.6, 0.025, 0.50))
    181    default_trees.append(tree("Lonicera",        1.0, 1.0, 7.0,  9.0, 2.0, 0.1, 0.6, 0.025, 0.70))
    182    default_trees.append(tree("Magnolia",        1.0, 1.0, 3.0,  5.0, 2.0, 0.1, 0.6, 0.025, 0.60))
    183    default_trees.append(tree("Malus",           1.0, 1.0, 4.5,  5.0, 2.0, 0.1, 0.6, 0.025, 0.30))
    184    default_trees.append(tree("Metasequoia",     5.0, 1.0, 4.5, 12.0, 2.0, 0.1, 0.6, 0.025, 0.50))
    185    default_trees.append(tree("Morus",           1.0, 1.0, 7.5, 11.5, 2.0, 0.1, 0.6, 0.025, 1.00))
    186    default_trees.append(tree("Ostrya",          1.0, 1.0, 2.0,  6.0, 2.0, 0.1, 0.6, 0.025, 1.00))
    187    default_trees.append(tree("Parrotia",        1.0, 1.0, 7.0,  7.0, 2.0, 0.1, 0.6, 0.025, 0.30))
    188    default_trees.append(tree("Paulownia",       1.0, 1.0, 4.0,  8.0, 2.0, 0.1, 0.6, 0.025, 0.40))
    189    default_trees.append(tree("Phellodendron",   1.0, 1.0,13.5, 13.5, 2.0, 0.1, 0.6, 0.025, 0.50))
    190    default_trees.append(tree("Picea",           3.0, 1.0, 3.0, 13.0, 2.0, 0.1, 0.6, 0.025, 0.90))
    191    default_trees.append(tree("Pinus",           3.0, 1.0, 6.0, 16.0, 2.0, 0.1, 0.6, 0.025, 0.80))
    192    default_trees.append(tree("Platanus",        1.0, 1.0,10.0, 14.5, 2.0, 0.1, 0.6, 0.025, 1.10))
    193    default_trees.append(tree("Populus",         1.0, 1.0, 9.0, 20.0, 2.0, 0.1, 0.6, 0.025, 1.40))
    194    default_trees.append(tree("Prunus",          1.0, 1.0, 5.0,  7.0, 2.0, 0.1, 0.6, 0.025, 1.60))
    195    default_trees.append(tree("Pseudotsuga",     3.0, 1.0, 6.0, 17.5, 2.0, 0.1, 0.6, 0.025, 0.70))
    196    default_trees.append(tree("Ptelea",          1.0, 1.0, 5.0,  4.0, 2.0, 0.1, 0.6, 0.025, 1.10))
    197    default_trees.append(tree("Pterocaria",      1.0, 1.0,10.0, 12.0, 2.0, 0.1, 0.6, 0.025, 0.50))
    198    default_trees.append(tree("Pterocarya",      1.0, 1.0,11.5, 14.5, 2.0, 0.1, 0.6, 0.025, 1.60))
    199    default_trees.append(tree("Pyrus",           3.0, 1.0, 3.0,  6.0, 2.0, 0.1, 0.6, 0.025, 1.80))
     135   default_trees.append(tree("Default",         1.0, 1.0, 4.0, 12.0, 3.0, 0.8, 0.6, 0.025, 0.35))
     136   default_trees.append(tree("Abies",           3.0, 1.0, 4.0, 12.0, 3.0, 0.8, 0.6, 0.025, 0.80))
     137   default_trees.append(tree("Acer",            1.0, 1.0, 7.0, 12.0, 3.0, 0.8, 0.6, 0.025, 0.80))
     138   default_trees.append(tree("Aesculus",        1.0, 1.0, 7.0, 12.0, 3.0, 0.8, 0.6, 0.025, 1.00))
     139   default_trees.append(tree("Ailanthus",       1.0, 1.0, 8.5, 13.5, 3.0, 0.8, 0.6, 0.025, 1.30))
     140   default_trees.append(tree("Alnus",           3.0, 1.0, 6.0, 16.0, 3.0, 0.8, 0.6, 0.025, 1.20))
     141   default_trees.append(tree("Amelanchier",     1.0, 1.0, 3.0,  4.0, 3.0, 0.8, 0.6, 0.025, 1.20))
     142   default_trees.append(tree("Betula",          1.0, 1.0, 6.0, 14.0, 3.0, 0.8, 0.6, 0.025, 0.30))
     143   default_trees.append(tree("Buxus",           1.0, 1.0, 4.0,  4.0, 3.0, 0.8, 0.6, 0.025, 0.90))
     144   default_trees.append(tree("Calocedrus",      3.0, 1.0, 5.0, 10.0, 3.0, 0.8, 0.6, 0.025, 0.50))
     145   default_trees.append(tree("Caragana",        1.0, 1.0, 3.5,  6.0, 3.0, 0.8, 0.6, 0.025, 0.90))
     146   default_trees.append(tree("Carpinus",        1.0, 1.0, 6.0, 10.0, 3.0, 0.8, 0.6, 0.025, 0.70))
     147   default_trees.append(tree("Carya",           1.0, 1.0, 5.0, 17.0, 3.0, 0.8, 0.6, 0.025, 0.80))
     148   default_trees.append(tree("Castanea",        1.0, 1.0, 4.5,  7.0, 3.0, 0.8, 0.6, 0.025, 0.80))
     149   default_trees.append(tree("Catalpa",         1.0, 1.0, 5.5,  6.5, 3.0, 0.8, 0.6, 0.025, 0.70))
     150   default_trees.append(tree("Cedrus",          1.0, 1.0, 8.0, 13.0, 3.0, 0.8, 0.6, 0.025, 0.80))
     151   default_trees.append(tree("Celtis",          1.0, 1.0, 6.0,  9.0, 3.0, 0.8, 0.6, 0.025, 0.80))
     152   default_trees.append(tree("Cercidiphyllum",  1.0, 1.0, 3.0,  6.5, 3.0, 0.8, 0.6, 0.025, 0.80))
     153   default_trees.append(tree("Cercis",          1.0, 1.0, 2.5,  7.5, 3.0, 0.8, 0.6, 0.025, 0.90))
     154   default_trees.append(tree("Chamaecyparis",   5.0, 1.0, 3.5,  9.0, 3.0, 0.8, 0.6, 0.025, 0.70))
     155   default_trees.append(tree("Cladrastis",      1.0, 1.0, 5.0, 10.0, 3.0, 0.8, 0.6, 0.025, 0.80))
     156   default_trees.append(tree("Cornus",          1.0, 1.0, 4.5,  6.5, 3.0, 0.8, 0.6, 0.025, 1.20))
     157   default_trees.append(tree("Corylus",         1.0, 1.0, 5.0,  9.0, 3.0, 0.8, 0.6, 0.025, 0.40))
     158   default_trees.append(tree("Cotinus",         1.0, 1.0, 4.0,  4.0, 3.0, 0.8, 0.6, 0.025, 0.70))
     159   default_trees.append(tree("Crataegus",       3.0, 1.0, 3.5,  6.0, 3.0, 0.8, 0.6, 0.025, 1.40))
     160   default_trees.append(tree("Cryptomeria",     3.0, 1.0, 5.0, 10.0, 3.0, 0.8, 0.6, 0.025, 0.50))
     161   default_trees.append(tree("Cupressocyparis", 3.0, 1.0, 3.0,  8.0, 3.0, 0.8, 0.6, 0.025, 0.40))
     162   default_trees.append(tree("Cupressus",       3.0, 1.0, 5.0,  7.0, 3.0, 0.8, 0.6, 0.025, 0.40))
     163   default_trees.append(tree("Cydonia",         1.0, 1.0, 2.0,  3.0, 3.0, 0.8, 0.6, 0.025, 0.90))
     164   default_trees.append(tree("Davidia",         1.0, 1.0,10.0, 14.0, 3.0, 0.8, 0.6, 0.025, 0.40))
     165   default_trees.append(tree("Elaeagnus",       1.0, 1.0, 6.5,  6.0, 3.0, 0.8, 0.6, 0.025, 1.20))
     166   default_trees.append(tree("Euodia",          1.0, 1.0, 4.5,  6.0, 3.0, 0.8, 0.6, 0.025, 0.90))
     167   default_trees.append(tree("Euonymus",        1.0, 1.0, 4.5,  6.0, 3.0, 0.8, 0.6, 0.025, 0.60))
     168   default_trees.append(tree("Fagus",           1.0, 1.0,10.0, 12.5, 3.0, 0.8, 0.6, 0.025, 0.50))
     169   default_trees.append(tree("Fraxinus",        1.0, 1.0, 5.5, 10.5, 3.0, 0.8, 0.6, 0.025, 1.60))
     170   default_trees.append(tree("Ginkgo",          3.0, 1.0, 4.0,  8.5, 3.0, 0.8, 0.6, 0.025, 0.80))
     171   default_trees.append(tree("Gleditsia",       1.0, 1.0, 6.5, 10.5, 3.0, 0.8, 0.6, 0.025, 0.60))
     172   default_trees.append(tree("Gymnocladus",     1.0, 1.0, 5.5, 10.0, 3.0, 0.8, 0.6, 0.025, 0.80))
     173   default_trees.append(tree("Hippophae",       1.0, 1.0, 9.5,  8.5, 3.0, 0.8, 0.6, 0.025, 0.80))
     174   default_trees.append(tree("Ilex",            1.0, 1.0, 4.0,  7.5, 3.0, 0.8, 0.6, 0.025, 0.80))
     175   default_trees.append(tree("Juglans",         1.0, 1.0, 7.0,  9.0, 3.0, 0.8, 0.6, 0.025, 0.50))
     176   default_trees.append(tree("Juniperus",       5.0, 1.0, 3.0,  7.0, 3.0, 0.8, 0.6, 0.025, 0.90))
     177   default_trees.append(tree("Koelreuteria",    1.0, 1.0, 3.5,  5.5, 3.0, 0.8, 0.6, 0.025, 0.50))
     178   default_trees.append(tree("Laburnum",        1.0, 1.0, 3.0,  6.0, 3.0, 0.8, 0.6, 0.025, 0.60))
     179   default_trees.append(tree("Larix",           3.0, 1.0, 7.0, 16.5, 3.0, 0.8, 0.6, 0.025, 0.60))
     180   default_trees.append(tree("Ligustrum",       1.0, 1.0, 3.0,  6.0, 3.0, 0.8, 0.6, 0.025, 1.10))
     181   default_trees.append(tree("Liquidambar",     3.0, 1.0, 3.0,  7.0, 3.0, 0.8, 0.6, 0.025, 0.30))
     182   default_trees.append(tree("Liriodendron",    3.0, 1.0, 4.5,  9.5, 3.0, 0.8, 0.6, 0.025, 0.50))
     183   default_trees.append(tree("Lonicera",        1.0, 1.0, 7.0,  9.0, 3.0, 0.8, 0.6, 0.025, 0.70))
     184   default_trees.append(tree("Magnolia",        1.0, 1.0, 3.0,  5.0, 3.0, 0.8, 0.6, 0.025, 0.60))
     185   default_trees.append(tree("Malus",           1.0, 1.0, 4.5,  5.0, 3.0, 0.8, 0.6, 0.025, 0.30))
     186   default_trees.append(tree("Metasequoia",     5.0, 1.0, 4.5, 12.0, 3.0, 0.8, 0.6, 0.025, 0.50))
     187   default_trees.append(tree("Morus",           1.0, 1.0, 7.5, 11.5, 3.0, 0.8, 0.6, 0.025, 1.00))
     188   default_trees.append(tree("Ostrya",          1.0, 1.0, 2.0,  6.0, 3.0, 0.8, 0.6, 0.025, 1.00))
     189   default_trees.append(tree("Parrotia",        1.0, 1.0, 7.0,  7.0, 3.0, 0.8, 0.6, 0.025, 0.30))
     190   default_trees.append(tree("Paulownia",       1.0, 1.0, 4.0,  8.0, 3.0, 0.8, 0.6, 0.025, 0.40))
     191   default_trees.append(tree("Phellodendron",   1.0, 1.0,13.5, 13.5, 3.0, 0.8, 0.6, 0.025, 0.50))
     192   default_trees.append(tree("Picea",           3.0, 1.0, 3.0, 13.0, 3.0, 0.8, 0.6, 0.025, 0.90))
     193   default_trees.append(tree("Pinus",           3.0, 1.0, 6.0, 16.0, 3.0, 0.8, 0.6, 0.025, 0.80))
     194   default_trees.append(tree("Platanus",        1.0, 1.0,10.0, 14.5, 3.0, 0.8, 0.6, 0.025, 1.10))
     195   default_trees.append(tree("Populus",         1.0, 1.0, 9.0, 20.0, 3.0, 0.8, 0.6, 0.025, 1.40))
     196   default_trees.append(tree("Prunus",          1.0, 1.0, 5.0,  7.0, 3.0, 0.8, 0.6, 0.025, 1.60))
     197   default_trees.append(tree("Pseudotsuga",     3.0, 1.0, 6.0, 17.5, 3.0, 0.8, 0.6, 0.025, 0.70))
     198   default_trees.append(tree("Ptelea",          1.0, 1.0, 5.0,  4.0, 3.0, 0.8, 0.6, 0.025, 1.10))
     199   default_trees.append(tree("Pterocaria",      1.0, 1.0,10.0, 12.0, 3.0, 0.8, 0.6, 0.025, 0.50))
     200   default_trees.append(tree("Pterocarya",      1.0, 1.0,11.5, 14.5, 3.0, 0.8, 0.6, 0.025, 1.60))
     201   default_trees.append(tree("Pyrus",           3.0, 1.0, 3.0,  6.0, 3.0, 0.8, 0.6, 0.025, 1.80))
    200202   default_trees.append(tree("Quercus",         1.0, 1.0, 8.0, 14.0, 3.1, 0.1, 0.6, 0.025, 0.40)) #
    201    default_trees.append(tree("Rhamnus",         1.0, 1.0, 4.5,  4.5, 2.0, 0.1, 0.6, 0.025, 1.30))
    202    default_trees.append(tree("Rhus",            1.0, 1.0, 7.0,  5.5, 2.0, 0.1, 0.6, 0.025, 0.50))
    203    default_trees.append(tree("Robinia",         1.0, 1.0, 4.5, 13.5, 2.0, 0.1, 0.6, 0.025, 0.50))
    204    default_trees.append(tree("Salix",           1.0, 1.0, 7.0, 14.0, 2.0, 0.1, 0.6, 0.025, 1.10))
    205    default_trees.append(tree("Sambucus",        1.0, 1.0, 8.0,  6.0, 2.0, 0.1, 0.6, 0.025, 1.40))
    206    default_trees.append(tree("Sasa",            1.0, 1.0,10.0, 25.0, 2.0, 0.1, 0.6, 0.025, 0.60))
    207    default_trees.append(tree("Sequoiadendron",  5.0, 1.0, 5.5, 10.5, 2.0, 0.1, 0.6, 0.025, 1.60))
    208    default_trees.append(tree("Sophora",         1.0, 1.0, 7.5, 10.0, 2.0, 0.1, 0.6, 0.025, 1.40))
    209    default_trees.append(tree("Sorbus",          1.0, 1.0, 4.0,  7.0, 2.0, 0.1, 0.6, 0.025, 1.10))
    210    default_trees.append(tree("Syringa",         1.0, 1.0, 4.5,  5.0, 2.0, 0.1, 0.6, 0.025, 0.60))
    211    default_trees.append(tree("Tamarix",         1.0, 1.0, 6.0,  7.0, 2.0, 0.1, 0.6, 0.025, 0.50))
    212    default_trees.append(tree("Taxodium",        5.0, 1.0, 6.0, 16.5, 2.0, 0.1, 0.6, 0.025, 0.60))
    213    default_trees.append(tree("Taxus",           2.0, 1.0, 5.0,  7.5, 2.0, 0.1, 0.6, 0.025, 1.50))
    214    default_trees.append(tree("Thuja",           3.0, 1.0, 3.5,  9.0, 2.0, 0.1, 0.6, 0.025, 0.70))
    215    default_trees.append(tree("Tilia",           3.0, 1.0, 7.0, 12.5, 2.0, 0.1, 0.6, 0.025, 0.70))
    216    default_trees.append(tree("Tsuga",           3.0, 1.0, 6.0, 10.5, 2.0, 0.1, 0.6, 0.025, 1.10))
    217    default_trees.append(tree("Ulmus",           1.0, 1.0, 7.5, 14.0, 2.0, 0.1, 0.6, 0.025, 0.80))
    218    default_trees.append(tree("Zelkova",         1.0, 1.0, 4.0,  5.5, 2.0, 0.1, 0.6, 0.025, 1.20))
    219    default_trees.append(tree("Zenobia",         1.0, 1.0, 5.0,  5.0, 2.0, 0.1, 0.6, 0.025, 0.40))
     203   default_trees.append(tree("Rhamnus",         1.0, 1.0, 4.5,  4.5, 3.0, 0.8, 0.6, 0.025, 1.30))
     204   default_trees.append(tree("Rhus",            1.0, 1.0, 7.0,  5.5, 3.0, 0.8, 0.6, 0.025, 0.50))
     205   default_trees.append(tree("Robinia",         1.0, 1.0, 4.5, 13.5, 3.0, 0.8, 0.6, 0.025, 0.50))
     206   default_trees.append(tree("Salix",           1.0, 1.0, 7.0, 14.0, 3.0, 0.8, 0.6, 0.025, 1.10))
     207   default_trees.append(tree("Sambucus",        1.0, 1.0, 8.0,  6.0, 3.0, 0.8, 0.6, 0.025, 1.40))
     208   default_trees.append(tree("Sasa",            1.0, 1.0,10.0, 25.0, 3.0, 0.8, 0.6, 0.025, 0.60))
     209   default_trees.append(tree("Sequoiadendron",  5.0, 1.0, 5.5, 10.5, 3.0, 0.8, 0.6, 0.025, 1.60))
     210   default_trees.append(tree("Sophora",         1.0, 1.0, 7.5, 10.0, 3.0, 0.8, 0.6, 0.025, 1.40))
     211   default_trees.append(tree("Sorbus",          1.0, 1.0, 4.0,  7.0, 3.0, 0.8, 0.6, 0.025, 1.10))
     212   default_trees.append(tree("Syringa",         1.0, 1.0, 4.5,  5.0, 3.0, 0.8, 0.6, 0.025, 0.60))
     213   default_trees.append(tree("Tamarix",         1.0, 1.0, 6.0,  7.0, 3.0, 0.8, 0.6, 0.025, 0.50))
     214   default_trees.append(tree("Taxodium",        5.0, 1.0, 6.0, 16.5, 3.0, 0.8, 0.6, 0.025, 0.60))
     215   default_trees.append(tree("Taxus",           2.0, 1.0, 5.0,  7.5, 3.0, 0.8, 0.6, 0.025, 1.50))
     216   default_trees.append(tree("Thuja",           3.0, 1.0, 3.5,  9.0, 3.0, 0.8, 0.6, 0.025, 0.70))
     217   default_trees.append(tree("Tilia",           3.0, 1.0, 7.0, 12.5, 3.0, 0.8, 0.6, 0.025, 0.70))
     218   default_trees.append(tree("Tsuga",           3.0, 1.0, 6.0, 10.5, 3.0, 0.8, 0.6, 0.025, 1.10))
     219   default_trees.append(tree("Ulmus",           1.0, 1.0, 7.5, 14.0, 3.0, 0.8, 0.6, 0.025, 0.80))
     220   default_trees.append(tree("Zelkova",         1.0, 1.0, 4.0,  5.5, 3.0, 0.8, 0.6, 0.025, 1.20))
     221   default_trees.append(tree("Zenobia",         1.0, 1.0, 5.0,  5.0, 3.0, 0.8, 0.6, 0.025, 0.40))
    220222   
    221223#  Define fill values
    … …  
    237239      if (season == "summer"):
    238240         tree_lai = default_trees[tree_type].lai_summer
     241
    239242      else:
    240243         tree_lai = default_trees[tree_type].lai_winter     
    … …  
    279282   
    280283   lad_max = tree_lai / (lad_max_part_1[0] + lad_max_part_2[0])
    281 
    282284   
    283285   
    … …  
    298300   z = np.arange(0,nz*dz,dz)
    299301   z[1:] = z[1:] - 0.5 * dz
    300 
    301302
    302303#  Define center of the tree position inside the local LAD domain
    … …  
    468469
    469470
    470 def process_patch(dz,patch_height,patch_lai,alpha,beta):
     471def process_patch(dz,patch_height,max_height_lad,patch_lai,alpha,beta):
    471472   
    472473#  Define fill values
    … …  
    481482   pch_index = np.where( (pch_index[:,:] == -1) ,0,pch_index[:,:])
    482483
    483    max_canopy_height = max(patch_height.flatten())
     484   max_canopy_height = max(max(patch_height.flatten()),max_height_lad)
    484485
    485486   z = np.arange(0,math.floor(max_canopy_height/dz)*dz+2*dz,dz)
    … …  
    500501          int_bpdf = 0.0
    501502          if ( (patch_height[j,i] != fillvalues["tree_data"]) & (patch_height[j,i] >= (0.5*dz)) ):
    502              for k in range(0,pch_index[j,i]):
     503             for k in range(1,pch_index[j,i]):
    503504                int_bpdf = int_bpdf + ( ( ( z[k] / patch_height[j,i] )**( alpha - 1 ) ) * ( ( 1.0 - ( z[k] / patch_height[j,i] ) )**(beta - 1 ) ) * ( dz / patch_height[j,i] ) )
    504505
    505              for k in range(0,pch_index[j,i]):
     506             for k in range(1,pch_index[j,i]):
    506507                pre_lad[k] =  patch_lai[j,i] * ( ( ( dz*k / patch_height[j,i] )**( alpha - 1.0 ) ) * ( ( 1.0 - ( dz*k / patch_height[j,i] ) )**(beta - 1.0 ) ) / int_bpdf ) / patch_height[j,i]
    507508
    … …  
    538539# dbh: default trunk diameter at breast height (1.4 m) (m)
    539540#
    540 #
    541541class tree:
    542542    def __init__(self, species, shape, ratio, diameter, height, lai_summer, lai_winter, lad_max_height, bad_scale, dbh):

  • palm/trunk/SCRIPTS/palm_csd_files/palm_csd_netcdf_interface.py

    r3629 r3668  
    2525# -----------------
    2626# $Id$
     27# Some improvements and new routines
     28#
     29# 3629 2018-12-13 12:18:54Z maronga
    2730# Some new routines
    2831#
    … …  
    3942
    4043from netCDF4 import Dataset
     44
     45def nc_read_from_file_1d_all(filename, varname):
     46
     47
     48   import numpy as np
     49   import sys
     50   
     51   try:
     52      f = open(filename)
     53      f.close()
     54#      print("Load: " + filename + ".")
     55   except FileNotFoundError:
     56      print("Error: " + filename + ". No such file. Aborting...")
     57      sys.exit(1)
     58   
     59   nc_file = Dataset(filename, "r+", format="NETCDF4")
     60   tmp_array = np.array(nc_file.variables[varname][:] , dtype=type(nc_file.variables[varname]))
     61
     62   return tmp_array
    4163
    4264def nc_read_from_file_2d(filename, varname, x0, x1, y0, y1):
    … …  
    174196   return 0
    175197
     198def nc_write_global_attribute(filename,attribute,value):
     199
     200   import datetime
     201   import os
     202
     203   print("Writing attribute " + attribute + " to file...")
     204
     205   f =  Dataset(filename, "a", format="NETCDF4")
     206   
     207   f.setncattr(attribute,value)
     208
     209   f.close()
     210
    176211
    177212def nc_write_dimension(filename,varname,array,datatype):

  • palm/trunk/SCRIPTS/palm_csd_files/palm_csd_tools.py

    r3567 r3668  
    2525# -----------------
    2626# $Id$
     27# Minor changes
     28#
     29# 3567 2018-11-27 13:59:21Z maronga
    2730# Initial revisions
    28 #
    29 #
    30 #
    3131#
    3232#
    … …  
    7777
    7878def interpolate_2d(array,x1,y1,x2,y2):
    79    
    80       tmp_int2d = interp2d(y1,x1,array,kind='linear')
    81       array_ip = tmp_int2d(y2.astype(float), x2.astype(float))
     79   tmp_int2d = interp2d(x1,y1,array,kind='linear')
     80   array_ip = tmp_int2d(x2.astype(float), y2.astype(float))
    8281
    83       return array_ip
     82   return array_ip
    8483   
    8584   
    … …  
    136135   return result
    137136
     137def check_consistency_3(array1,array2,array3,fill1,fill2,fill3):
     138
     139   tmp_array = np.where(array1 != fill1,1,0) + np.where(array2 != fill2,1,0) + np.where(array3 != fill3,-1,0)
     140   
     141   test = np.any(tmp_array.flatten() != 0)
     142   if test:
     143      print("soil_type array is not consistent!")
     144      print("max: " + str(max(tmp_array.flatten())) + ", min: " + str(min(tmp_array.flatten())))
     145   else:
     146      print("soil_type array is consistent!")     
     147   return tmp_array, test
     148
     149
     150
    138151
    139152def check_consistency_4(array1,array2,array3,array4,fill1,fill2,fill3,fill4):
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |