Changeset 3726 for palm/trunk/SCRIPTS


Ignore:
Timestamp:
Feb 7, 2019 6:22:49 PM (5 years ago)
Author:
maronga
Message:

bugfixes in palm_csd

Location:
palm/trunk/SCRIPTS
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SCRIPTS/palm_csd

    r3688 r3726  
    2525# -----------------
    2626# $Id$
     27# Removed some more bugs
     28#
     29# 3688 2019-01-22 10:44:20Z maronga
    2730# Some unspecified bugfixes
    2831#
     
    458461   filename.append(path_out + "/" + filename_out + "_" + domain_names[i])
    459462   if domain_parent[i] is not None:
    460       ii_parent.append(domain_names.index(domain_parent[i]))
     463      ii_parent.append(input_px.index(domain_px[domain_names.index(domain_parent[i])]))
    461464   else:
    462465      ii_parent.append(None)
    463466
    464467
    465    x_UTM = nc_read_from_file_2d(input_file_x[ii[i]], "Band1", domain_x0[i], domain_x0[i], domain_y0[i], domain_y0[i]) - 0.5 * domain_px[i]
    466    y_UTM = nc_read_from_file_2d(input_file_y[ii[i]], "Band1", domain_x0[i], domain_x0[i], domain_y0[i], domain_y0[i]) - 0.5 * domain_px[i]
    467    lat = nc_read_from_file_2d(input_file_lat[ii[i]], "Band1", domain_x0[i], domain_x0[i], domain_y0[i], domain_y0[i])
    468    lon = nc_read_from_file_2d(input_file_lon[ii[i]], "Band1", domain_x0[i], domain_x0[i], domain_y0[i], domain_y0[i]) 
     468   x_UTM = nc_read_from_file_2d(input_file_x[ii[i]], "Band1", domain_x0[i], domain_x0[i]+1, domain_y0[i], domain_y0[i]+1)
     469   y_UTM = nc_read_from_file_2d(input_file_y[ii[i]], "Band1", domain_x0[i], domain_x0[i]+1, domain_y0[i], domain_y0[i]+1)
     470   lat = nc_read_from_file_2d(input_file_lat[ii[i]], "Band1", domain_x0[i], domain_x0[i]+1, domain_y0[i], domain_y0[i]+1)
     471   lon = nc_read_from_file_2d(input_file_lon[ii[i]], "Band1", domain_x0[i], domain_x0[i]+1, domain_y0[i], domain_y0[i]+1) 
     472   
     473#  Calculate position of origin
     474   x_UTM_origin = float(x_UTM[0,0]) - 0.5 * (float(x_UTM[0,1]) - float(x_UTM[0,0]))
     475   y_UTM_origin = float(y_UTM[0,0]) - 0.5 * (float(y_UTM[1,0]) - float(y_UTM[0,0])) 
     476   x_origin = float(lon[0,0]) - 0.5 * (float(lon[0,1]) - float(lon[0,0]))
     477   y_origin = float(lat[0,0]) - 0.5 * (float(lat[1,0]) - float(lat[0,0]))   
    469478   
    470479#  Create NetCDF output file and set global attributes
    471480   nc_create_file(filename[i])
    472    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)
     481   nc_write_global_attributes(filename[i],x_UTM_origin,y_UTM_origin,y_origin,x_origin,"",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)
    473482
    474483   del x_UTM, y_UTM, lat, lon
     
    492501del zt_all[:]
    493502
    494 print( "Shift terrain height by -" + str(zt_min))
     503print( "Shift terrain heights by -" + str(zt_min))
    495504for i in range(0,ndomains):
    496505
     
    500509   y = nc_read_from_file_1d(input_file_y[ii[i]], "y", domain_y0[i], domain_y1[i])
    501510
     511 
    502512   zt = zt - zt_min
    503513   
     
    506516#  If necessary, interpolate parent domain terrain height on child domain grid and blend the two
    507517   if domain_ip[i]:
    508       tmp_x0 = int( domain_x0[i] * domain_px[i] / domain_px[ii_parent[i]] ) - 1
    509       tmp_y0 = int( domain_y0[i] * domain_px[i] / domain_px[ii_parent[i]] ) - 1
    510       tmp_x1 = int( domain_x1[i] * domain_px[i] / domain_px[ii_parent[i]] ) + 1
    511       tmp_y1 = int( domain_y1[i] * domain_px[i] / domain_px[ii_parent[i]] ) + 1
    512                
     518      parent_id = domain_names.index(domain_parent[i])
     519      tmp_x0 = int( domain_x0[i] * domain_px[i] / domain_px[parent_id] ) - 1
     520      tmp_y0 = int( domain_y0[i] * domain_px[i] / domain_px[parent_id] ) - 1
     521      tmp_x1 = int( domain_x1[i] * domain_px[i] / domain_px[parent_id] ) + 1
     522      tmp_y1 = int( domain_y1[i] * domain_px[i] / domain_px[parent_id] ) + 1
     523      print(tmp_x0)
     524      print(tmp_x1)     
     525      print(tmp_y0)
     526      print(tmp_y1) 
     527
    513528      tmp_x = nc_read_from_file_1d(input_file_x[ii_parent[i]], "x", tmp_x0, tmp_x1)
    514529      tmp_y = nc_read_from_file_1d(input_file_y[ii_parent[i]], "y", tmp_y0, tmp_y1)
     
    516531      zt_parent = nc_read_from_file_2d(input_file_zt[ii_parent[i]], 'Band1', tmp_x0, tmp_x1, tmp_y0, tmp_y1)
    517532
    518       print( "Shift terrain height by -" + str(zt_min))
    519533      zt_parent = zt_parent - zt_min
    520534
    521535#     Interpolate array and bring to PALM grid of child domain
    522 
    523536      zt_ip = interpolate_2d(zt_parent,tmp_x,tmp_y,x,y)
    524       zt_ip = bring_to_palm_grid(zt_ip,x,y,domain_dz[ii_parent[i]])
    525 
    526 #     Shift the child terrain height according to the parent mean terrain height       
    527       zt = zt - np.mean(zt)  + np.mean(zt_ip)
     537      zt_ip = bring_to_palm_grid(zt_ip,x,y,domain_dz[parent_id])
     538     
     539     
     540#     Shift the child terrain height according to the parent mean terrain height
     541      zt = zt - np.mean(zt) + np.mean(zt_ip)
    528542 
    529543 
     
    532546   
    533547#  Final step: add zt array to the global array   
     548
    534549   zt_all.append(zt)
    535550   del zt
     
    569584   x_UTM = nc_read_from_file_2d(input_file_x_UTM[ii[i]], "Band1", domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])
    570585   y_UTM = nc_read_from_file_2d(input_file_y_UTM[ii[i]], "Band1", domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) 
     586
    571587   
    572588   nc_write_to_file_2d(filename[i], 'E_UTM', x_UTM, datatypes["E_UTM"],'y','x',fillvalues["E_UTM"])
     
    664680      building_id = np.where(bridges_2d == fillvalues["bridges_2d"],building_id,bridges_id)
    665681     
    666      
    667682      if np.any(buildings_2d != fillvalues["buildings_2d"]):
    668683         buildings_3d, z = make_3d_from_2d(buildings_2d,x,y,domain_dz[i])
     
    676691         nc_write_attribute(filename[i], 'z', 'long_name', 'z') 
    677692         nc_write_attribute(filename[i], 'z', 'units', 'm')
     693         
     694         nc_overwrite_to_file_2d(filename[i], 'building_id', building_id) 
    678695 
    679696         nc_write_to_file_3d(filename[i], 'buildings_3d', buildings_3d, datatypes["buildings_3d"],'z','y','x',fillvalues["buildings_3d"])   
     
    839856   street_crossings = np.where((street_crossings < 1) & (street_crossings != fillvalues["street_crossings"]),defaultvalues["street_crossings"],street_crossings)
    840857   
    841    nc_write_to_file_2d(filename[i], 'street_crossings', street_crossings, datatypes["street_crossings"],'y','x',fillvalues["street_crossings"])
    842    nc_write_attribute(filename[i], 'street_crossings', 'long_name', 'street crossings')
    843    nc_write_attribute(filename[i], 'street_crossings', 'units', '')
    844    nc_write_attribute(filename[i], 'street_crossings', 'res_orig', domain_px[i]) 
    845    nc_write_attribute(filename[i], 'street_crossings', 'coordinates', 'E_UTM N_UTM lon lat')
    846    nc_write_attribute(filename[i], 'street_crossings', 'grid_mapping', 'E_UTM N_UTM lon lat')   
     858   nc_write_to_file_2d(filename[i], 'street_crossing', street_crossings, datatypes["street_crossings"],'y','x',fillvalues["street_crossings"])
     859   nc_write_attribute(filename[i], 'street_crossing', 'long_name', 'street crossings')
     860   nc_write_attribute(filename[i], 'street_crossing', 'units', '')
     861   nc_write_attribute(filename[i], 'street_crossing', 'res_orig', domain_px[i]) 
     862   nc_write_attribute(filename[i], 'street_crossing', 'coordinates', 'E_UTM N_UTM lon lat')
     863   nc_write_attribute(filename[i], 'street_crossing', 'grid_mapping', 'E_UTM N_UTM lon lat')   
    847864   del street_crossings
    848865
     
    10171034      vegetation_pars = nc_read_from_file_3d_all(filename[i], 'vegetation_pars')
    10181035      lai = vegetation_pars[1,:,:]
    1019 
    1020 #     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
     1036     
     1037
     1038#     Determine all pixels that do not already have an LAD but which are high vegetation to a dummy value of 1.0 and remove all other pixels
    10211039      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"])
    10221040
    1023 #     Now, assign either the default LAI for high vegetation or the value stored in the LAI array.
    1024       lai_high = np.where( (lai_high != fillvalues["tree_data"]) & (lai == fillvalues["tree_data"]), settings_lai_high_default, lai)
     1041#     Now, assign either the default LAI for high vegetation or keep 1.0 from the lai_high array.
     1042      lai_high = np.where( (lai_high != fillvalues["tree_data"]) & (lai == fillvalues["tree_data"]), settings_lai_high_default, lai_high)
     1043
     1044#     If lai values are available in the lai array, write them on the lai_high array
     1045      lai_high = np.where( (lai_high != fillvalues["tree_data"]) & (lai != fillvalues["tree_data"]), lai, lai_high)
    10251046
    10261047#     Define a patch height wherever it is missing, but where a high vegetation LAI was set
     
    10381059         
    10391060      if ( max(patch_height.flatten()) >= (2.0 * domain_dz[i]) ):
    1040          print("    start calculating LAD (this might take some time)")
     1061         print("    start calculating LAD (this might take some time)")       
     1062
     1063         
    10411064         lad_patch, patch_nz, status = process_patch(domain_dz[i],patch_height,max(zlad),lai_high,settings_lai_alpha,settings_lai_beta)
    10421065
  • palm/trunk/SCRIPTS/palm_csd_files/palm_csd_tools.py

    r3668 r3726  
    2525# -----------------
    2626# $Id$
     27# Index bound bugfix
     28#
     29# 3668 2019-01-14 12:49:24Z maronga
    2730# Minor changes
    2831#
     
    9093      for l in range(0,len(x)):
    9194         for m in range(0,len(y)):
    92             for k in range(0,len(k_tmp+1)):
     95            for k in range(1,len(k_tmp+1)):
    9396               if k_tmp[k] > array[m,l]:
    94                   array[m,l] = k_tmp[k]-0.5*dz
     97                  array[m,l] = k_tmp[k-1]
    9598                  break
    9699
     
    105108      array_3d = np.ones((len(k_tmp),len(y),len(x)))
    106109     
    107       for l in range(0,len(x)-1):
    108          for m in range(0,len(y)-1):
    109             for k in range(0,len(k_tmp)-1):
     110      for l in range(0,len(x)):
     111         for m in range(0,len(y)):
     112            for k in range(0,len(k_tmp)):
    110113               if k_tmp[k] > array_2d[m,l]:
    111114                  array_3d[k,m,l] = 0
Note: See TracChangeset for help on using the changeset viewer.