Changeset 4763 for palm/trunk/SCRIPTS/palm_cvd

Timestamp:

Oct 30, 2020 12:06:47 PM (4 years ago)

Author:

suehring

Message:

Revision of setting-up virtual measurements. virtual_measurement_mod: Simplified input of coordinates. All coordinate and height arrays are now 1D, independent on featuretype. This allows easier usage also for other campaign data sets independent of UC2; Avoid type conversion from 64 to 32 bit when output the data; Increase dimension size of sample-variable string (sometimes more than 100 variables are listed for heavy measured locations); Remove quantities that cannot be sampled; | Further, revision of the script palm_cvd to convert measurement coordinates from UC2 data standard towards a PALM readable format: Convert trajectory and timeseriesProfile coordinates into 1-D coordinates equivalent to timeseries coordiates. This simplifies processing in PALM and makes the virtual-measurement module also applicable to other campaigns; Check automatically for data organization (stored in subdirectories or not).

File:

• palm/trunk/SCRIPTS/palm_cvd (modified) (9 diffs)

palm/trunk/SCRIPTS/palm_cvd

@@ -26 +26 @@
 # -----------------
 # $Id$
+# - Check automatically for data organization (stored in subdirectories or not)
+# - Convert trajectory and timeseriesProfile coordinates into 1-D coordinates
+#   equivalent to timeseries coordiates. This simplifies processing in PALM and
+#   makes the virtual-measurement module also applicable to other campaigns.
+#
+# 4758 2020-10-26 13:03:52Z suehring
 # In order to do not omit observations that are on the same site but have different
 # coordinates or feature-types, process all files rather than only one and omit 
@@ -54 +60 @@
 import os
 import numpy as np
+import time
 
 
@@ -187 +194 @@
 name_featuretype   = "featureType"
 name_ts            = "timeSeries"
+name_tspr          = "timeSeriesProfile"
 name_traj          = "trajectory"
 name_ntime         = "ntime"
@@ -196 +204 @@
 name_eutm          = "E_UTM"
 name_nutm          = "N_UTM"
-name_hao           = "height_above_origin"
+name_hao           = "height"
 name_station_h     = "station_h"
 name_z             = "z"
@@ -271 +279 @@
 if ( input_from_observations ):
 
-   # Run loop over all subdirectories, detect the files and extract a list of sites.
+   # Run loop over all listed input data. Depending on the data set, this could be
+   # a list of files or a list of subdirectories.
    # This is done to reduce the number of virtual measurements in the model. Each
    # virtual measurement has an overhead and consumes memory.
    sites = []
    input_files = []
+   input_files_orig = []
    for dirname in list_input_data:
        data_file = data_path + dirname
 
-       # Directory may contain various file versions.
-       # Take the one with highest cycle number.
-       highest_cycle_nr = 0
-       for filename in os.listdir(data_file):
-          start_seq = len( filename ) - 6
-          end_seq   = len( filename ) - 3
-          if int( filename[start_seq:end_seq] ) > highest_cycle_nr:
-             highest_cycle_nr = int(filename[start_seq:end_seq])
-             latest_file      = filename
+       if ( os.path.isdir(data_file) == True ):
+          # Directory may contain various file versions.
+          # Take the one with highest cycle number.
+          highest_cycle_nr = 0
+          for filename in os.listdir(data_file):
+             start_seq = len( filename ) - 6
+             end_seq   = len( filename ) - 3
+             if int( filename[start_seq:end_seq] ) > highest_cycle_nr:
+                highest_cycle_nr = int(filename[start_seq:end_seq])
+                latest_file      = filename
+          input_file = data_file + "/" + latest_file
+          input_file_orig = latest_file
+       else:
+          input_file = data_file
+          input_file_orig = dirname
 
        # Open the NetCDF file
-       input_file = data_file + "/" + latest_file
        ncfile_in = Dataset( input_file, "r", format="NETCDF4", encoding='ascii')
        input_files.append(input_file)
-
-   # Define a nested list of default variables that shall be measured. Based on this list,
-   # the final number of measured variables is determined.
-   measured_variables_all_sites = [ ['u', 'v', 'w', 'theta', 'hus'] for var in range(0, len(input_files))]
-
-   # Run loop over all subdirectories that contain observational data
-   for counter, file in enumerate(input_files):
-
-      # Open the NetCDF input file
-      input_file = input_files[counter]
+       input_files_orig.append(input_file_orig)
+
+#  Gather all files according to their feature type and all sites for the respective feature type
+   files_traj  = []
+   files_ts    = []
+   files_tspr  = []
+   sites_traj  = []
+   sites_ts    = []
+   sites_tspr  = []
+   for input_file in input_files:
       ncfile_in = Dataset( input_file, "r", format="NETCDF4", encoding='ascii' )
 
-      # Determine index for the treated site
-      num_vmeas = input_files.index( input_file ) + 1
-      print( counter, num_vmeas )
-
-      # Read global attributes and write them immediately into the output file
       for att in ncfile_in.ncattrs():
          if ( att == name_featuretype ):
             feature = ncfile_in.getncattr(att)
-         if ( att == name_datacontent ):
-            content = ncfile_in.getncattr(att)
          if ( att == name_site ):
             site = ncfile_in.getncattr(att)
 
-         if ( att in atts_float ):
-            ncfile_out.setncattr( att + str(num_vmeas), np.double(ncfile_in.getncattr(att)) )
-         else:
-            ncfile_out.setncattr( att + str(num_vmeas), ncfile_in.getncattr(att) )
-
-      #timeSeries
-      if ( feature == name_ts ):
-         ntime = len( ncfile_in.dimensions[name_ntime]   )
-         nstat = len( ncfile_in.dimensions[name_station] )
-         ncfile_out.createDimension( name_ntime   + str(num_vmeas), ntime )
-         ncfile_out.createDimension( name_station + str(num_vmeas), nstat )
-
-      #trajectory
-      elif ( feature == name_traj ):
-         ntime = len( ncfile_in.dimensions[name_ntime]   )
-         ntraj = len( ncfile_in.dimensions[name_traj_dim] )
-         ncfile_out.createDimension( name_ntime    + str(num_vmeas), ntime )
-         ncfile_out.createDimension( name_traj_dim + str(num_vmeas), ntraj )
-
-      #timeseriesProfile
+      if ( feature == name_traj ):
+         files_traj.append(input_file)
+      elif ( feature == name_ts ):
+         files_ts.append(input_file)
       else:
-         ntime = len( ncfile_in.dimensions[name_ntime]   )
-         nstat = len( ncfile_in.dimensions[name_station] )
-         nz    = len( ncfile_in.dimensions[name_nz]      )
-         ncfile_out.createDimension( name_ntime   + str(num_vmeas), ntime )
-         ncfile_out.createDimension( name_station + str(num_vmeas), nstat )
-         ncfile_out.createDimension( name_nz      + str(num_vmeas), nz    )
-
-      for var in ncfile_in.variables.keys():
-         if ( var in dims_out ):
-            # Create a variable and write it to file after it is read. In order to
-            # avoid fill values in the dimensions, these are converted to zero
-            # before written to file. Depending on the featureType of the measurement,
-            # the array shape is different. For more informations, please see
-            # [UC]2 data standard.
-            # Timeseries
-            if ( feature == name_ts  ):
-               temp_ts = ncfile_out.createVariable( var + str(num_vmeas), float, \
-                                                    name_station + str(num_vmeas))
-               temp_ts[:] = np.nan_to_num( ncfile_in.variables[var][:] )
-
-            # Trajectories
-            elif ( feature == name_traj ):
-               # @note: If there are files where station_h is present although featureType is
-               #        trajectory, station_h must not be read.
-               if var != name_station_h:
-                  temp_traj = ncfile_out.createVariable( var + str(num_vmeas), float, \
-                                                         ( name_traj_dim + str(num_vmeas), \
-                                                           name_ntime + str(num_vmeas) ) )
-                  temp_traj[:,:] = np.nan_to_num( ncfile_in.variables[var][:,:] )
-
-            # TimeseriesProfiles
-            else:
-               if ( var == 'z' ):
-                  temp_pr = ncfile_out.createVariable( var + str(num_vmeas), float, \
-                                                      ( name_station + str(num_vmeas), \
-                                                        name_nz + str(num_vmeas) ) )
-                  temp_pr[:] = np.nan_to_num( ncfile_in.variables[var][:,0,:] )
-               else:
-                  temp_pr = ncfile_out.createVariable( var + str(num_vmeas), float, \
-                                                       name_station + str(num_vmeas))
-                  temp_pr[:] = np.nan_to_num( ncfile_in.variables[var][:] )
-
-      # Search for variables to be measured. In case the variable isn't already defined,
-      # append the variable to the list.
-      for var in ncfile_in.variables.keys():
-         if ( var not in non_measurable_vars  and  \
-              var not in vars_default         and  \
-              var not in measured_variables_all_sites[input_files.index( input_file )] ):
-
-            measured_variables_all_sites[input_files.index( input_file )].append(var)
-
-      # Close the NetCDF input file
+         files_tspr.append(input_file)
+
+      if ( feature == name_traj  and  site not in sites_traj ):
+         sites_traj.append(site)
+      if ( feature == name_ts  and  site not in sites_ts ):
+         sites_ts.append(site)
+      if ( feature == name_tspr  and  site not in sites_tspr ):
+         sites_tspr.append(site)
+
       ncfile_in.close()
 
-   # After variables are gathered and dimensions / attributes are already written to file,
-   # the list of measured variables is written to file.
-   for site in input_files:
-
-      num_vmeas = input_files.index( site ) + 1
-
-      ncfile_out.createDimension( "nvar"+ str(num_vmeas), \
-                                  len( measured_variables_all_sites[input_files.index( site )] ) )
-
-      measured = ncfile_out.createVariable( 'measured_variables' + str(num_vmeas), 'S1', \
-                                            ("nvar" + str(num_vmeas), "string_len")) # must be NC_CHAR
-
-      for counter, meas in enumerate( measured_variables_all_sites[input_files.index( site )] ):
-         measured[counter] = stringtochar( np.array( meas,"S%s"%(max_string_len) ) )
-
-      # Check if any of the measured variables is a soil variable. Set flag accordingly.
+   for input_file in files_traj:
+      print( "traj", input_file )
+   for site in sites_traj:
+      print( "traj", site )
+
+   for site in sites_ts:
+      print( "ts", site )
+   for site in sites_tspr:
+      print( "tspr", site )
+      
+   for file in files_tspr:
+      print( "tspr", file )
+   counter_id = 1
+   for site_traj in sites_traj:
+      # For the given site already define the featureTpye and site
+      ncfile_out.setncattr( name_featuretype + str(counter_id), name_traj )
+      ncfile_out.setncattr( name_site + str(counter_id), site_traj )
+
+      # Define the number of coordinates for the site
+      num_coords = 0
+
+      e_utm_traj = np.array([])
+      n_utm_traj = np.array([])
+      h_traj     = np.array([])
+      measured_variables = ['u', 'v', 'w', 'theta', 'hus']
+      for input_file in files_traj:
+         ncfile_in = Dataset( input_file, "r", format="NETCDF4", encoding='ascii' )
+         for att in ncfile_in.ncattrs():
+            if ( att == name_site ):
+               site = ncfile_in.getncattr(att)
+
+         if ( site == site_traj ):
+
+            for att in ncfile_in.ncattrs():
+               if ( att == name_orig_x ):
+                  orig_x = ncfile_in.getncattr(att)
+               if ( att == name_orig_y ):
+                  orig_y = ncfile_in.getncattr(att)
+               if ( att == name_orig_z ):
+                  orig_z = ncfile_in.getncattr(att)
+
+            ntime = len( ncfile_in.dimensions[name_ntime]    )
+            ntraj = len( ncfile_in.dimensions[name_traj_dim] )
+
+            num_coords += ntime * ntraj
+#           Gather UTM and height coordinates and merge them into one array. Further, gather
+#           the variables that shall be sampled. Coordinates are checked to for NaN values and
+#           are tranformed to arithmetric numbers. Further, 2D input array is transformed into
+#           a 1D array.
+            for var in ncfile_in.variables.keys():
+               if ( var in dims_out  and  var == name_eutm ):
+                  e_utm_traj = np.append(e_utm_traj, np.nan_to_num( ncfile_in.variables[var][:,:] ).flatten())
+                  #e_utm_traj.append( np.nan_to_num( ncfile_in.variables[var][:,:] ).flatten() )
+               if ( var in dims_out  and  var == name_nutm ):
+                  n_utm_traj = np.append(n_utm_traj, np.nan_to_num( ncfile_in.variables[var][:,:] ).flatten())
+               if ( var in dims_out  and  var == name_hao ):
+                  h_traj = np.append(h_traj, np.nan_to_num( ncfile_in.variables[var][:,:] ).flatten())
+
+               if ( var not in non_measurable_vars  and  \
+                    var not in vars_default         and  \
+                    var not in measured_variables ):
+                  measured_variables.append(var)
+
+         ncfile_in.close()
+
+#     After all files for the current site are processed, write the origin-coordinates for x,y,z
+      ncfile_out.setncattr( name_orig_x + str(counter_id), orig_x )
+      ncfile_out.setncattr( name_orig_y + str(counter_id), orig_y )
+      ncfile_out.setncattr( name_orig_z + str(counter_id), orig_z )
+#     Create the dimensions
+      ncfile_out.createDimension( name_station + str(counter_id), num_coords )
+
+      temp_traj = ncfile_out.createVariable( name_eutm + str(counter_id), float, name_station + str(counter_id) )
+      temp_traj[:] = e_utm_traj
+
+      temp_traj = ncfile_out.createVariable( name_nutm + str(counter_id), float, name_station + str(counter_id) )
+      temp_traj[:] = n_utm_traj
+
+      temp_traj = ncfile_out.createVariable( name_hao  + str(counter_id), float, name_station + str(counter_id) )
+      temp_traj[:] = h_traj
+
+#     Check if any of the measured variables is a soil variable. Set flag accordingly.
       soil = False
-      for var in measured_variables_all_sites[input_files.index( site )]:
+      for var in measured_variables:
          if ( var in soil_vars ):
             soil = True
-
-      # Write soil flag
-      ncfile_out.setncattr( name_soil_sampling + str( num_vmeas), np.int8(soil) )
-
-
-
-    # Store the number of observational sites
-   num_sites += len( input_files )
+#     Write soil flag
+      ncfile_out.setncattr( name_soil_sampling + str( counter_id), np.int8(soil) )
+
+#     Create dimension for sample-variable string
+      ncfile_out.createDimension( "nvar"+ str(counter_id), len( measured_variables ) )
+
+      measured_var = ncfile_out.createVariable( 'measured_variables' + str(counter_id), 'S1', \
+                                                ("nvar" + str(counter_id), "string_len") ) # must be NC_CHAR
+
+      # Write the variables to the file
+      for counter, meas in enumerate( measured_variables ):
+         measured_var[counter] = stringtochar( np.array( meas,"S%s"%(max_string_len) ) )
+
+#     Increment the counter
+      counter_id += 1
+   
+   
+   for site_tspr in sites_tspr:
+      # For the given site already define the featureTpye and site
+      ncfile_out.setncattr( name_featuretype + str(counter_id), name_tspr )
+      ncfile_out.setncattr( name_site + str(counter_id), site_tspr )
+
+      # Define the number of coordinates for the site
+      num_coords = 0
+      e_utm_tspr     = np.array([])
+      n_utm_tspr     = np.array([])
+      station_h_tspr = np.array([])
+      z_tspr         = np.array([])
+
+      measured_variables = ['u', 'v', 'w', 'theta', 'hus']
+      for input_file in files_tspr:
+         ncfile_in = Dataset( input_file, "r", format="NETCDF4", encoding='ascii' )
+         for att in ncfile_in.ncattrs():
+            if ( att == name_site ):
+               site = ncfile_in.getncattr(att)
+
+         if ( site == site_tspr ):
+            for att in ncfile_in.ncattrs():
+               if ( att == name_orig_x ):
+                  orig_x = ncfile_in.getncattr(att)
+               if ( att == name_orig_y ):
+                  orig_y = ncfile_in.getncattr(att)
+               if ( att == name_orig_z ):
+                  orig_z = ncfile_in.getncattr(att)
+
+            nstation = len( ncfile_in.dimensions[name_station] )
+            ntime    = len( ncfile_in.dimensions[name_ntime]   )
+            nz       = len( ncfile_in.dimensions[name_nz]   )
+
+            num_coords += nstation * ntime * nz
+#           Gather UTM and height coordinates and merge them into one array. Further, gather
+#           the variables that shall be sampled. Coordinates are checked to for NaN values and
+#           are tranformed to arithmetric numbers. Further, 2D input array is transformed into
+#           a 1D array.
+            for var in ncfile_in.variables.keys():
+               tspr_tmp1 = np.zeros((nstation))
+               tspr_tmp2 = np.zeros((ntime*nz))
+               if ( var in dims_out  and  var == name_eutm ):
+                  tspr_tmp1 = np.nan_to_num( ncfile_in.variables[var][:] )
+                  for ns in range(0,int(nstation)):
+                     tspr_tmp2[:] = tspr_tmp1[ns]
+                     e_utm_tspr = np.append(e_utm_tspr, tspr_tmp2)
+               if ( var in dims_out  and  var == name_nutm ):
+                  tspr_tmp1 = np.nan_to_num( ncfile_in.variables[var][:] )
+                  for ns in range(0,int(nstation)):
+                     tspr_tmp2[:] = tspr_tmp1[ns]
+                     n_utm_tspr = np.append(n_utm_tspr, tspr_tmp2)
+               if ( var in dims_out  and  var == name_z ):
+                  z_tspr_tmp = np.nan_to_num( ncfile_in.variables[var][:,:,:] )
+                  z_tspr = np.append(z_tspr, np.concatenate( z_tspr_tmp ))
+               if ( var in dims_out  and  var == name_station_h ):
+                  tspr_tmp1 = np.nan_to_num( ncfile_in.variables[var][:] )
+                  for ns in range(0,int(nstation)):
+                     tspr_tmp2[:] = tspr_tmp1[ns]
+                     station_h_tspr = np.append(station_h_tspr, tspr_tmp2)
+
+               if ( var not in non_measurable_vars  and  \
+                    var not in vars_default         and  \
+                    var not in measured_variables ):
+                  measured_variables.append(var)
+
+         ncfile_in.close()
+
+#     After all files for the current site are processed, write the origin-coordinates for x,y,z
+      ncfile_out.setncattr( name_orig_x + str(counter_id), orig_x )
+      ncfile_out.setncattr( name_orig_y + str(counter_id), orig_y )
+      ncfile_out.setncattr( name_orig_z + str(counter_id), orig_z )
+#     Create the dimensions
+      ncfile_out.createDimension( name_station + str(counter_id), num_coords )
+
+      temp_tspr = ncfile_out.createVariable( name_eutm      + str(counter_id), float, name_station + str(counter_id) )
+      temp_tspr[:] = e_utm_tspr
+
+      temp_tspr = ncfile_out.createVariable( name_nutm      + str(counter_id), float, name_station + str(counter_id) )
+      temp_tspr[:] = n_utm_tspr
+
+      temp_tspr = ncfile_out.createVariable( name_z         + str(counter_id), float, name_station + str(counter_id) )
+      temp_tspr[:] = z_tspr
+
+      temp_tspr = ncfile_out.createVariable( name_station_h + str(counter_id), float, name_station + str(counter_id) )
+      temp_tspr[:] = station_h_tspr
+
+#     Check if any of the measured variables is a soil variable. Set flag accordingly.
+      soil = False
+      for var in measured_variables:
+         if ( var in soil_vars ):
+            soil = True
+#     Write soil flag
+      ncfile_out.setncattr( name_soil_sampling + str( counter_id), np.int8(soil) )
+
+#     Create dimension for sample-variable string
+      ncfile_out.createDimension( "nvar"+ str(counter_id), len( measured_variables ) )
+
+      measured_var = ncfile_out.createVariable( 'measured_variables' + str(counter_id), 'S1', \
+                                                ("nvar" + str(counter_id), "string_len") ) # must be NC_CHAR
+
+      # Write the variables to the file
+      for counter, meas in enumerate( measured_variables ):
+         measured_var[counter] = stringtochar( np.array( meas,"S%s"%(max_string_len) ) )
+
+#     Increment the counter
+      counter_id += 1
+   
+   
+   for site_ts in sites_ts:
+      # For the given site already define the featureTpye and site
+      ncfile_out.setncattr( name_featuretype + str(counter_id), name_ts )
+      ncfile_out.setncattr( name_site + str(counter_id), site_ts )
+
+      # Define the number of coordinates for the site
+      num_coords = 0
+      e_utm_ts     = np.array([])
+      n_utm_ts     = np.array([])
+      station_h_ts = np.array([])
+      z_ts         = np.array([])
+      
+      measured_variables = ['u', 'v', 'w', 'theta', 'hus']
+      for input_file in files_ts:
+         ncfile_in = Dataset( input_file, "r", format="NETCDF4", encoding='ascii' )
+         for att in ncfile_in.ncattrs():
+            if ( att == name_site ):
+               site = ncfile_in.getncattr(att)
+
+         if ( site == site_ts ):
+
+            for att in ncfile_in.ncattrs():
+               if ( att == name_orig_x ):
+                  orig_x = ncfile_in.getncattr(att)
+               if ( att == name_orig_y ):
+                  orig_y = ncfile_in.getncattr(att)
+               if ( att == name_orig_z ):
+                  orig_z = ncfile_in.getncattr(att)
+
+            nstation = len( ncfile_in.dimensions[name_station]    )
+            num_coords += nstation
+#           Gather UTM and height coordinates and merge them into one array. Further, gather
+#           the variables that shall be sampled. Coordinates are checked to for NaN values and
+#           are tranformed to arithmetric numbers.
+            for var in ncfile_in.variables.keys():
+               if ( var in dims_out  and  var == name_eutm ):
+                  e_utm_ts = np.append(e_utm_ts, np.nan_to_num( ncfile_in.variables[var][:] ))
+               if ( var in dims_out  and  var == name_nutm ):
+                  n_utm_ts = np.append(n_utm_ts, np.nan_to_num( ncfile_in.variables[var][:] ))
+               if ( var in dims_out  and  var == name_z ):
+                  z_ts = np.append(z_ts, np.nan_to_num( ncfile_in.variables[var][:] ))
+               if ( var in dims_out  and  var == name_station_h ):
+                  station_h_ts = np.append(station_h_ts, np.nan_to_num( ncfile_in.variables[var][:] ))
+
+               if ( var not in non_measurable_vars  and  \
+                    var not in vars_default         and  \
+                    var not in measured_variables ):
+                  measured_variables.append(var)
+
+         ncfile_in.close()
+
+#     After all files for the current site are processed, write the origin-coordinates for x,y,z
+      ncfile_out.setncattr( name_orig_x + str(counter_id), orig_x )
+      ncfile_out.setncattr( name_orig_y + str(counter_id), orig_y )
+      ncfile_out.setncattr( name_orig_z + str(counter_id), orig_z )
+#     Create the dimensions
+      ncfile_out.createDimension( name_station + str(counter_id), num_coords )
+
+      temp_ts = ncfile_out.createVariable( name_eutm      + str(counter_id), float, name_station + str(counter_id) )
+      temp_ts[:] = e_utm_ts
+
+      temp_ts = ncfile_out.createVariable( name_nutm      + str(counter_id), float, name_station + str(counter_id) )
+      temp_ts[:] = n_utm_ts
+
+      temp_ts = ncfile_out.createVariable( name_z         + str(counter_id), float, name_station + str(counter_id) )
+      temp_ts[:] = z_ts
+
+      temp_ts = ncfile_out.createVariable( name_station_h + str(counter_id), float, name_station + str(counter_id) )
+      temp_ts[:] = station_h_ts
+
+#     Check if any of the measured variables is a soil variable. Set flag accordingly.
+      soil = False
+      for var in measured_variables:
+         if ( var in soil_vars ):
+            soil = True
+#     Write soil flag
+      ncfile_out.setncattr( name_soil_sampling + str( counter_id), np.int8(soil) )
+
+#     Create dimension for sample-variable string
+      ncfile_out.createDimension( "nvar"+ str(counter_id), len( measured_variables ) )
+
+      measured_var = ncfile_out.createVariable( 'measured_variables' + str(counter_id), 'S1', \
+                                                ("nvar" + str(counter_id), "string_len") ) # must be NC_CHAR
+
+      # Write the variables to the file
+      for counter, meas in enumerate( measured_variables ):
+         measured_var[counter] = stringtochar( np.array( meas,"S%s"%(max_string_len) ) )
+
+#     Increment the counter
+      counter_id += 1
+
+# Store the number of observational sites
+   num_sites = len( sites_traj ) + len( sites_ts ) + len( sites_tspr )
 
 
@@ -426 +650 @@
 # are possible.
 if ( custom_coordinates ):
-
-   count_site = num_sites + 1
+   num_sites = counter_id - 1
    for coord in coordinates:
       # Define mandatory attributes
-      ncfile_out.setncattr( name_featuretype + str(count_site),  \
+      ncfile_out.setncattr( name_featuretype + str(counter_id),  \
                             name_ts )
-      ncfile_out.setncattr( name_site        + str(count_site),  \
-                            "custom"         + str(count_site - num_sites) )
-      ncfile_out.setncattr( name_orig_x      + str(count_site),  \
+      ncfile_out.setncattr( name_site        + str(counter_id),  \
+                            "custom"         + str(counter_id - num_sites) )
+      ncfile_out.setncattr( name_orig_x      + str(counter_id),  \
                             coord[0] )
-      ncfile_out.setncattr( name_orig_y      + str(count_site),  \
+      ncfile_out.setncattr( name_orig_y      + str(counter_id),  \
                             coord[1] )
-      ncfile_out.setncattr( name_orig_z      + str(count_site),  \
+      ncfile_out.setncattr( name_orig_z      + str(counter_id),  \
                             0.0 )
 
@@ -444 +667 @@
       ntime = 1
       nstat = 1
-      ncfile_out.createDimension( name_ntime   + str(count_site), ntime )
-      ncfile_out.createDimension( name_station + str(count_site), nstat )
+      ncfile_out.createDimension( name_ntime   + str(counter_id), ntime )
+      ncfile_out.createDimension( name_station + str(counter_id), nstat )
 
       # Define coordinate variables
-      temp_ts = ncfile_out.createVariable( name_eutm      + str(count_site), \
+      temp_ts = ncfile_out.createVariable( name_eutm      + str(counter_id), \
                                            float,                            \
-                                           name_station   + str(count_site) )
+                                           name_station   + str(counter_id) )
       temp_ts[:] = np.array( coord[0] )
 
-      temp_ts = ncfile_out.createVariable( name_nutm      + str(count_site), \
+      temp_ts = ncfile_out.createVariable( name_nutm      + str(counter_id), \
                                            float,                            \
-                                           name_station   + str(count_site) )
+                                           name_station   + str(counter_id) )
       temp_ts[:] = np.array( coord[1] )
 
-      temp_ts = ncfile_out.createVariable( name_z         + str(count_site), \
+      temp_ts = ncfile_out.createVariable( name_z         + str(counter_id), \
                                            float,                            \
-                                           name_station   + str(count_site) )
+                                           name_station   + str(counter_id) )
       temp_ts[:] = np.array( coord[2] )
 
-      temp_ts = ncfile_out.createVariable( name_station_h + str(count_site), \
+      temp_ts = ncfile_out.createVariable( name_station_h + str(counter_id), \
                                            float,                            \
-                                           name_station   + str(count_site) )
+                                           name_station   + str(counter_id) )
       temp_ts[:] = np.array( 0.0 )
 
 
-      count_site += 1
+      counter_id += 1
 
    # Reset counter variable
-   count_site = num_sites + 1
+   counter_id = num_sites + 1
 
    # check if variables are prescribed. If so, prepare final output string
@@ -491 +714 @@
                 measured_variables.append(var)
 
-         ncfile_out.createDimension( "nvar"+ str(count_site), \
+         ncfile_out.createDimension( "nvar"+ str(counter_id), \
                                      len( measured_variables ) )
 
-         measured_var = ncfile_out.createVariable( 'measured_variables' + str(count_site), 'S1', \
-                                                  ("nvar" + str(count_site), "string_len") ) # must be NC_CHAR
+         measured_var = ncfile_out.createVariable( 'measured_variables' + str(counter_id), 'S1', \
+                                                  ("nvar" + str(counter_id), "string_len") ) # must be NC_CHAR
 
          # Write the variables to the file
@@ -507 +730 @@
 
          # Write soil flag
-         ncfile_out.setncattr( name_soil_sampling + str( count_site), np.int8(soil) )
+         ncfile_out.setncattr( name_soil_sampling + str( counter_id), np.int8(soil) )
 
          # Increment counter variable
-         count_site += 1
+         counter_id += 1
 
          del ( measured_variables[:] )
 
    # Add the number of customized sites.
-   num_sites += int( number_positions )
+   num_sites = counter_id - 1