Changeset 3726 for palm/trunk/SCRIPTS

Timestamp:

Feb 7, 2019 6:22:49 PM (6 years ago)

Author:

maronga

Message:

bugfixes in palm_csd

Location:

palm/trunk/SCRIPTS

Files:

• palm_csd (modified) (13 diffs)
• palm_csd_files/palm_csd_tools.py (modified) (3 diffs)

palm/trunk/SCRIPTS/palm_csd

@@ -25 +25 @@
 # -----------------
 # $Id$
+# Removed some more bugs
+# 
+# 3688 2019-01-22 10:44:20Z maronga
 # Some unspecified bugfixes
 # 
@@ -458 +461 @@
    filename.append(path_out + "/" + filename_out + "_" + domain_names[i])
    if domain_parent[i] is not None:
-      ii_parent.append(domain_names.index(domain_parent[i]))
+      ii_parent.append(input_px.index(domain_px[domain_names.index(domain_parent[i])]))
    else:
       ii_parent.append(None)
 
 
-   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]
-   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]
-   lat = nc_read_from_file_2d(input_file_lat[ii[i]], "Band1", domain_x0[i], domain_x0[i], domain_y0[i], domain_y0[i])
-   lon = nc_read_from_file_2d(input_file_lon[ii[i]], "Band1", domain_x0[i], domain_x0[i], domain_y0[i], domain_y0[i])  
+   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)
+   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)
+   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)
+   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)  
+   
+#  Calculate position of origin
+   x_UTM_origin = float(x_UTM[0,0]) - 0.5 * (float(x_UTM[0,1]) - float(x_UTM[0,0]))
+   y_UTM_origin = float(y_UTM[0,0]) - 0.5 * (float(y_UTM[1,0]) - float(y_UTM[0,0]))  
+   x_origin = float(lon[0,0]) - 0.5 * (float(lon[0,1]) - float(lon[0,0]))
+   y_origin = float(lat[0,0]) - 0.5 * (float(lat[1,0]) - float(lat[0,0]))   
    
 #  Create NetCDF output file and set global attributes
    nc_create_file(filename[i])
-   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)
+   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)
 
    del x_UTM, y_UTM, lat, lon
@@ -492 +501 @@
 del zt_all[:]
 
-print( "Shift terrain height by -" + str(zt_min))
+print( "Shift terrain heights by -" + str(zt_min))
 for i in range(0,ndomains):
 
@@ -500 +509 @@
    y = nc_read_from_file_1d(input_file_y[ii[i]], "y", domain_y0[i], domain_y1[i])
 
+ 
    zt = zt - zt_min
    
@@ -506 +516 @@
 #  If necessary, interpolate parent domain terrain height on child domain grid and blend the two
    if domain_ip[i]:
-      tmp_x0 = int( domain_x0[i] * domain_px[i] / domain_px[ii_parent[i]] ) - 1
-      tmp_y0 = int( domain_y0[i] * domain_px[i] / domain_px[ii_parent[i]] ) - 1
-      tmp_x1 = int( domain_x1[i] * domain_px[i] / domain_px[ii_parent[i]] ) + 1
-      tmp_y1 = int( domain_y1[i] * domain_px[i] / domain_px[ii_parent[i]] ) + 1
-                
+      parent_id = domain_names.index(domain_parent[i])
+      tmp_x0 = int( domain_x0[i] * domain_px[i] / domain_px[parent_id] ) - 1
+      tmp_y0 = int( domain_y0[i] * domain_px[i] / domain_px[parent_id] ) - 1
+      tmp_x1 = int( domain_x1[i] * domain_px[i] / domain_px[parent_id] ) + 1
+      tmp_y1 = int( domain_y1[i] * domain_px[i] / domain_px[parent_id] ) + 1
+      print(tmp_x0)
+      print(tmp_x1)      
+      print(tmp_y0)
+      print(tmp_y1)  
+
       tmp_x = nc_read_from_file_1d(input_file_x[ii_parent[i]], "x", tmp_x0, tmp_x1)
       tmp_y = nc_read_from_file_1d(input_file_y[ii_parent[i]], "y", tmp_y0, tmp_y1)
@@ -516 +531 @@
       zt_parent = nc_read_from_file_2d(input_file_zt[ii_parent[i]], 'Band1', tmp_x0, tmp_x1, tmp_y0, tmp_y1)
 
-      print( "Shift terrain height by -" + str(zt_min))
       zt_parent = zt_parent - zt_min
 
 #     Interpolate array and bring to PALM grid of child domain
-
       zt_ip = interpolate_2d(zt_parent,tmp_x,tmp_y,x,y)
-      zt_ip = bring_to_palm_grid(zt_ip,x,y,domain_dz[ii_parent[i]])
-
-#     Shift the child terrain height according to the parent mean terrain height       
-      zt = zt - np.mean(zt)  + np.mean(zt_ip)
+      zt_ip = bring_to_palm_grid(zt_ip,x,y,domain_dz[parent_id])
+      
+      
+#     Shift the child terrain height according to the parent mean terrain height 
+      zt = zt - np.mean(zt) + np.mean(zt_ip)
   
   
@@ -532 +546 @@
    
 #  Final step: add zt array to the global array   
+
    zt_all.append(zt)
    del zt
@@ -569 +584 @@
    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])
    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])  
+
    
    nc_write_to_file_2d(filename[i], 'E_UTM', x_UTM, datatypes["E_UTM"],'y','x',fillvalues["E_UTM"])
@@ -664 +680 @@
       building_id = np.where(bridges_2d == fillvalues["bridges_2d"],building_id,bridges_id)
       
-      
       if np.any(buildings_2d != fillvalues["buildings_2d"]):
          buildings_3d, z = make_3d_from_2d(buildings_2d,x,y,domain_dz[i])
@@ -676 +691 @@
          nc_write_attribute(filename[i], 'z', 'long_name', 'z')  
          nc_write_attribute(filename[i], 'z', 'units', 'm')
+         
+         nc_overwrite_to_file_2d(filename[i], 'building_id', building_id)  
  
          nc_write_to_file_3d(filename[i], 'buildings_3d', buildings_3d, datatypes["buildings_3d"],'z','y','x',fillvalues["buildings_3d"])    
@@ -839 +856 @@
    street_crossings = np.where((street_crossings < 1) & (street_crossings != fillvalues["street_crossings"]),defaultvalues["street_crossings"],street_crossings)
    
-   nc_write_to_file_2d(filename[i], 'street_crossings', street_crossings, datatypes["street_crossings"],'y','x',fillvalues["street_crossings"]) 
-   nc_write_attribute(filename[i], 'street_crossings', 'long_name', 'street crossings')
-   nc_write_attribute(filename[i], 'street_crossings', 'units', '')
-   nc_write_attribute(filename[i], 'street_crossings', 'res_orig', domain_px[i])  
-   nc_write_attribute(filename[i], 'street_crossings', 'coordinates', 'E_UTM N_UTM lon lat') 
-   nc_write_attribute(filename[i], 'street_crossings', 'grid_mapping', 'E_UTM N_UTM lon lat')   
+   nc_write_to_file_2d(filename[i], 'street_crossing', street_crossings, datatypes["street_crossings"],'y','x',fillvalues["street_crossings"]) 
+   nc_write_attribute(filename[i], 'street_crossing', 'long_name', 'street crossings')
+   nc_write_attribute(filename[i], 'street_crossing', 'units', '')
+   nc_write_attribute(filename[i], 'street_crossing', 'res_orig', domain_px[i])  
+   nc_write_attribute(filename[i], 'street_crossing', 'coordinates', 'E_UTM N_UTM lon lat') 
+   nc_write_attribute(filename[i], 'street_crossing', 'grid_mapping', 'E_UTM N_UTM lon lat')   
    del street_crossings
 
@@ -1017 +1034 @@
       vegetation_pars = nc_read_from_file_3d_all(filename[i], 'vegetation_pars') 
       lai = vegetation_pars[1,:,:]
-
-#     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
+     
+
+#     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
       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"])
 
-#     Now, assign either the default LAI for high vegetation or the value stored in the LAI array.
-      lai_high = np.where( (lai_high != fillvalues["tree_data"]) & (lai == fillvalues["tree_data"]), settings_lai_high_default, lai)
+#     Now, assign either the default LAI for high vegetation or keep 1.0 from the lai_high array.
+      lai_high = np.where( (lai_high != fillvalues["tree_data"]) & (lai == fillvalues["tree_data"]), settings_lai_high_default, lai_high)
+
+#     If lai values are available in the lai array, write them on the lai_high array
+      lai_high = np.where( (lai_high != fillvalues["tree_data"]) & (lai != fillvalues["tree_data"]), lai, lai_high)
 
 #     Define a patch height wherever it is missing, but where a high vegetation LAI was set
@@ -1038 +1059 @@
          
       if ( max(patch_height.flatten()) >= (2.0 * domain_dz[i]) ):
-         print("    start calculating LAD (this might take some time)") 
+         print("    start calculating LAD (this might take some time)")       
+
+         
          lad_patch, patch_nz, status = process_patch(domain_dz[i],patch_height,max(zlad),lai_high,settings_lai_alpha,settings_lai_beta)
 

palm/trunk/SCRIPTS/palm_csd_files/palm_csd_tools.py

@@ -25 +25 @@
 # -----------------
 # $Id$
+# Index bound bugfix
+# 
+# 3668 2019-01-14 12:49:24Z maronga
 # Minor changes
 # 
@@ -90 +93 @@
       for l in range(0,len(x)):
          for m in range(0,len(y)):
-            for k in range(0,len(k_tmp+1)):
+            for k in range(1,len(k_tmp+1)):
                if k_tmp[k] > array[m,l]:
-                  array[m,l] = k_tmp[k]-0.5*dz
+                  array[m,l] = k_tmp[k-1]
                   break
 
@@ -105 +108 @@
       array_3d = np.ones((len(k_tmp),len(y),len(x)))
       
-      for l in range(0,len(x)-1):
-         for m in range(0,len(y)-1):
-            for k in range(0,len(k_tmp)-1):
+      for l in range(0,len(x)):
+         for m in range(0,len(y)):
+            for k in range(0,len(k_tmp)):
                if k_tmp[k] > array_2d[m,l]:
                   array_3d[k,m,l] = 0