Merge pull request #249 from xylar/land_ice_mesh_tools_to_python3

Land ice mesh tools to python 3

Make a few minor formatting changes to the land ice mesh tools so they work in python 3 as well as python 2.

I also include mark_horns_for_culling.py because I think this script might be used primarily by the land-ice model.

Author: xylar

landice/mesh_tools_li/create_landice_grid_from_generic_MPAS_grid.py

@@ -3,6 +3,9 @@
 # I've only tested it with a periodic_hex grid, but it should work with any MPAS grid.
 # Currently variable attributes are not copied (and periodic_hex does not assign any, so this is ok).  If variable attributes are added to periodic_hex, this script should be modified to copy them (looping over dir(var), skipping over variable function names "assignValue", "getValue", "typecode").
 
+from __future__ import absolute_import, division, print_function, \
+    unicode_literals
+
 import sys, numpy
 from netCDF4 import Dataset
 from optparse import OptionParser
@@ -11,7 +14,7 @@
 
 sphere_radius = 6.37122e6 # earth radius, if needed
 
-print "** Gathering information.  (Invoke with --help for more details. All arguments are optional)"
+print("** Gathering information.  (Invoke with --help for more details. All arguments are optional)")
 parser = OptionParser()
 parser.add_option("-i", "--in", dest="fileinName", help="input filename.  Defaults to 'grid.nc'", metavar="FILENAME")
 parser.add_option("-o", "--out", dest="fileoutName", help="output filename.  Defaults to 'landice_grid.nc'", metavar="FILENAME")
@@ -26,14 +29,14 @@
 options, args = parser.parse_args()
 
 if not options.fileinName:
-    print "No input filename specified, so using 'grid.nc'."
+    print("No input filename specified, so using 'grid.nc'.")
     options.fileinName = 'grid.nc'
 else:
-    print "Input file is:", options.fileinName
+    print("Input file is: {}".format(options.fileinName))
 if not options.fileoutName:
-    print "No output filename specified, so using 'landice_grid.nc'."
+    print("No output filename specified, so using 'landice_grid.nc'.")
     options.fileoutName = 'landice_grid.nc'
-print '' # make a space in stdout before further output
+print('') # make a space in stdout before further output
 
 # Get the input file
 filein = Dataset(options.fileinName,'r')
@@ -47,27 +50,27 @@
 # ============================================
 # Do this first as doing it last is slow for big files since adding
 # attributes forces the contents to get reorganized.
-print "---- Copying global attributes from input file to output file ----"
+print("---- Copying global attributes from input file to output file ----")
 for name in filein.ncattrs():
   # sphere radius needs to be set to that of the earth if on a sphere
   if name == 'sphere_radius' and getattr(filein, 'on_a_sphere') == "YES             ":
     setattr(fileout, 'sphere_radius', sphere_radius)
-    print 'Set global attribute   sphere_radius = ', str(sphere_radius)
+    print('Set global attribute   sphere_radius = {}'.format(sphere_radius))
   elif name =='history':
     # Update history attribute of netCDF file
     newhist = '\n'.join([getattr(filein, 'history'), ' '.join(sys.argv[:]) ] )
     setattr(fileout, 'history', newhist )
   else:
     # Otherwise simply copy the attr
     setattr(fileout, name, getattr(filein, name) )
-    print 'Copied global attribute  ', name, '=', getattr(filein, name)
+    print('Copied global attribute  {} = {}'.format(name, getattr(filein, name)))
 
 # Update history attribute of netCDF file if we didn't above
 if not hasattr(fileout, 'history'):
    setattr(fileout, 'history', sys.argv[:] )
 
 fileout.sync()
-print '' # make a space in stdout before further output
+print('') # make a space in stdout before further output
 
 
 # ============================================
@@ -78,7 +81,7 @@
 #       It may be better to list them explicitly as I do for the grid variables,
 #       but this way ensures they all get included and is easier.
 # Note: The UNLIMITED time dimension will return a dimension value of None with Scientific.IO.  This is what is supposed to happen.  See below for how to deal with assigning values to a variable with a unlimited dimension.  Special handling is needed with the netCDF module.
-print "---- Copying dimensions from input file to output file ----"
+print("---- Copying dimensions from input file to output file ----")
 for dim in filein.dimensions.keys():
     if dim == 'nTracers':
         pass  # Do nothing - we don't want this dimension
@@ -91,12 +94,12 @@
         if options.levels is None:
           # If nVertLevels is in the input file, and a value for it was not
           # specified on the command line, then use the value from the file (do nothing here)
-          print "Using nVertLevels from the intput file:", len(filein.dimensions[dim])
+          print("Using nVertLevels from the intput file: {}".format(len(filein.dimensions[dim])))
           dimvalue = len(filein.dimensions[dim])
         else:
           # if nVertLevels is in the input file, but a value WAS specified
           # on the command line, then use the command line value
-          print "Using nVertLevels specified on the command line:", int(options.levels)
+          print("Using nVertLevels specified on the command line: {}".format(int(options.levels)))
           dimvalue = int(options.levels)
       else:
          dimvalue = len(filein.dimensions[dim])
@@ -105,31 +108,31 @@
 # it has not been added to the output file yet.  Treat those here.
 if 'nVertLevels' not in fileout.dimensions:
    if options.levels is None:
-       print "nVertLevels not in input file and not specified.  Using default value of 10."
+       print("nVertLevels not in input file and not specified.  Using default value of 10.")
        fileout.createDimension('nVertLevels', 10)
    else:
-       print "Using nVertLevels specified on the command line:", int(options.levels)
+       print("Using nVertLevels specified on the command line: {}".format(int(options.levels)))
        fileout.createDimension('nVertLevels', int(options.levels))
 # Also create the nVertInterfaces dimension, even if none of the variables require it.
 fileout.createDimension('nVertInterfaces', len(fileout.dimensions['nVertLevels']) + 1)  # nVertInterfaces = nVertLevels + 1
-print 'Added new dimension nVertInterfaces to output file with value of ' + str(len(fileout.dimensions['nVertInterfaces'])) + '.'
+print('Added new dimension nVertInterfaces to output file with value of {}.'.format(len(fileout.dimensions['nVertInterfaces'])))
 
 fileout.sync()
-print 'Finished creating dimensions in output file.\n' # include an extra blank line here
+print('Finished creating dimensions in output file.\n') # include an extra blank line here
 
 # ============================================
 # Copy over all of the required grid variables to the new file
 # ============================================
-print "Beginning to copy mesh variables to output file."
+print("Beginning to copy mesh variables to output file.")
 vars2copy = ['latCell', 'lonCell', 'xCell', 'yCell', 'zCell', 'indexToCellID', 'latEdge', 'lonEdge', 'xEdge', 'yEdge', 'zEdge', 'indexToEdgeID', 'latVertex', 'lonVertex', 'xVertex', 'yVertex', 'zVertex', 'indexToVertexID', 'cellsOnEdge', 'nEdgesOnCell', 'nEdgesOnEdge', 'edgesOnCell', 'edgesOnEdge', 'weightsOnEdge', 'dvEdge', 'dcEdge', 'angleEdge', 'areaCell', 'areaTriangle', 'cellsOnCell', 'verticesOnCell', 'verticesOnEdge', 'edgesOnVertex', 'cellsOnVertex', 'kiteAreasOnVertex']
 # Add these optional fields if they exist in the input file
 for optionalVar in ['meshDensity', 'gridSpacing', 'cellQuality', 'triangleQuality', 'triangleAngleQuality', 'obtuseTriangle']:
    if optionalVar in filein.variables:
       vars2copy.append(optionalVar)
 
 for varname in vars2copy:
-   print "-",
-print "|"
+   print("-"),
+print("|")
 for varname in vars2copy:
    thevar = filein.variables[varname]
    datatype = thevar.dtype
@@ -146,8 +149,8 @@
    del newVar, thevar
    sys.stdout.write("* "); sys.stdout.flush()
 fileout.sync()
-print "|"
-print "Finished copying mesh variables to output file.\n"
+print("|")
+print("Finished copying mesh variables to output file.\n")
 
 # ============================================
 # Create the land ice variables (all the shallow water vars in the input file can be ignored)
@@ -170,7 +173,7 @@
       layerInterfaces[k] = 4.0/3.0 * (1.0 - ((k+1.0-1.0)/(nInterfaces-1.0) + 1.0)**-2)
    for k in range(nVertLevels):
       layerThicknessFractionsData[k] = layerInterfaces[k+1] - layerInterfaces[k]
-   print "Setting layerThicknessFractions to:", layerThicknessFractionsData
+   print("Setting layerThicknessFractions to: {}".format(layerThicknessFractionsData))
 else:
    sys.exit('Unknown method for vertical spacing method (--vert): '+options.vertMethod)
 
@@ -192,17 +195,17 @@
 newvar[:] = numpy.zeros(newvar.shape)
 newvar = fileout.createVariable('floatingBasalMassBal', datatype, ('Time', 'nCells'))
 newvar[:] = numpy.zeros(newvar.shape)
-print 'Added default variables: thickness, temperature, bedTopography, sfcMassBal, floatingBasalMassBal'
+print('Added default variables: thickness, temperature, bedTopography, sfcMassBal, floatingBasalMassBal')
 
 if options.beta:
    newvar = fileout.createVariable('beta', datatype, ('Time', 'nCells'))
    newvar[:] = 1.0e8  # Give a default beta that won't have much sliding.
-   print 'Added optional variable: beta'
+   print('Added optional variable: beta')
 
 if options.effecpress:
    newvar = fileout.createVariable('effectivePressure', datatype, ('Time', 'nCells'))
    newvar[:] = 1.0e8  # Give a default effective pressure that won't have much sliding.
-   print 'Added optional variable: effectivePressure'
+   print('Added optional variable: effectivePressure')
 
 if options.dirichlet:
    newvar = fileout.createVariable('dirichletVelocityMask', datatypeInt, ('Time', 'nCells', 'nVertInterfaces'))
@@ -211,7 +214,7 @@
    newvar[:] = 0.0
    newvar = fileout.createVariable('uReconstructY', datatype, ('Time', 'nCells', 'nVertInterfaces',))
    newvar[:] = 0.0
-   print 'Added optional dirichlet variables: dirichletVelocityMask, uReconstructX, uReconstructY'
+   print('Added optional dirichlet variables: dirichletVelocityMask, uReconstructX, uReconstructY')
 
 if options.thermal:
    newvar = fileout.createVariable('temperature', datatype, ('Time', 'nCells', 'nVertLevels'))
@@ -220,7 +223,7 @@
    newvar[:] = 273.15 # Give default value for temperate ice
    newvar = fileout.createVariable('basalHeatFlux', datatype, ('Time', 'nCells'))
    newvar[:] = 0.0 # Default to none (W/m2)
-   print 'Added optional thermal variables: temperature, surfaceAirTemperature, basalHeatFlux'
+   print('Added optional thermal variables: temperature, surfaceAirTemperature, basalHeatFlux')
 
 if options.hydro:
    newvar = fileout.createVariable('waterThickness', datatype, ('Time', 'nCells'))
@@ -237,7 +240,7 @@
    newvar[:] = 0.0
    newvar = fileout.createVariable('waterFluxMask', 'i', ('Time', 'nEdges'))
    newvar[:] = 0.0
-   print 'Added optional hydro variables: waterThickness, tillWaterThickness, meltInput, frictionAngle, waterPressure, waterFluxMask'
+   print('Added optional hydro variables: waterThickness, tillWaterThickness, meltInput, frictionAngle, waterPressure, waterFluxMask')
 
 if options.obs:
    newvar = fileout.createVariable('observedSurfaceVelocityX', datatype, ('Time', 'nCells'))
@@ -252,7 +255,7 @@
    newvar[:] = 0.0
    newvar = fileout.createVariable('thicknessUncertainty', datatype, ('Time', 'nCells'))
    newvar[:] = 0.0
-   print 'Added optional velocity optimization variables: observedSurfaceVelocityX, observedSurfaceVelocityY, observedSurfaceVelocityUncertainty, observedThicknessTendency, observedThicknessTendencyUncertainty, thicknessUncertainty'
+   print('Added optional velocity optimization variables: observedSurfaceVelocityX, observedSurfaceVelocityY, observedSurfaceVelocityUncertainty, observedThicknessTendency, observedThicknessTendencyUncertainty, thicknessUncertainty')
 
 # Update history attribute of netCDF file
 thiscommand = datetime.now().strftime("%a %b %d %H:%M:%S %Y") + ": " + " ".join(sys.argv[:])
@@ -262,10 +265,10 @@
    newhist = thiscommand
 setattr(fileout, 'history', newhist )
 
-print "Completed creating land ice variables in new file. Now syncing to file."
+print("Completed creating land ice variables in new file. Now syncing to file.")
 fileout.sync()
 
 filein.close()
 fileout.close()
 
-print '\n** Successfully created ' + options.fileoutName + '.**'
+print('\n** Successfully created {}.**'.format(options.fileoutName))

landice/mesh_tools_li/define_cullMask.py

@@ -2,14 +2,17 @@
 # Script for adding a field named cullMask to an MPAS land ice grid for use with the MpasCellCuller tool that actually culls the unwanted cells.
 # Matt Hoffman, February 28, 2013
 
+from __future__ import absolute_import, division, print_function, \
+    unicode_literals
+
 import sys
 import numpy as np
 from optparse import OptionParser
 from netCDF4 import Dataset as NetCDFFile
 from datetime import datetime
 
 
-print "** Gathering information."
+print("** Gathering information.")
 parser = OptionParser()
 parser.add_option("-f", "--file", dest="file", help="grid file to modify; default: landice_grid.nc", metavar="FILE")
 parser.add_option("-m", "--method", dest="method", help="method to use for marking cells to cull.  Supported methods: 'noIce', 'numCells', 'distance', 'radius', 'edgeFraction'", metavar="METHOD")
@@ -19,7 +22,7 @@
 options, args = parser.parse_args()
 
 if not options.file:
-	print "No grid filename provided. Using landice_grid.nc."
+        print("No grid filename provided. Using landice_grid.nc.")
         options.file = "landice_grid.nc"
 
 if not options.method:
@@ -46,18 +49,18 @@
 thicknessMissing = True
 try:
   thickness = f.variables['thickness'][0,:]
-  print 'Using thickness field at time 0'
+  print('Using thickness field at time 0')
   thicknessMissing = False
 except:
-  print "The field 'thickness' is not available.  Some culling methods will not work."
+  print("The field 'thickness' is not available.  Some culling methods will not work.")
 
 
 # =====  Various methods for defining the mask ====
 
 # =========
 #  only keep cells with ice
 if maskmethod == 'noIce':
-  print "Method: remove cells without ice"
+  print("Method: remove cells without ice")
   if thicknessMissing:
      sys.exit("Unable to perform 'numCells' method because thickness field was missing.")
 
@@ -66,13 +69,13 @@
 # =========
 #  add a buffer of X cells around the ice
 elif maskmethod == 'numCells':
-  print "Method: remove cells beyond a certain number of cells from existing ice"
+  print("Method: remove cells beyond a certain number of cells from existing ice")
 
   if thicknessMissing:
      sys.exit("Unable to perform 'numCells' method because thickness field was missing.")
 
   buffersize=int(options.numCells)  # number of cells to expand
-  print "Using a buffer of {} cells".format(buffersize)
+  print("Using a buffer of {} cells".format(buffersize))
 
   keepCellMask = np.copy(cullCell[:])
   keepCellMask[:] = 0
@@ -81,17 +84,17 @@
 
   # mark the cells with ice first
   keepCellMask[thickness > 0.0] = 1
-  print 'Num of cells with ice:', sum(keepCellMask)
+  print('Num of cells with ice: {}'.format(sum(keepCellMask)))
 
   for i in range(buffersize):
-    print 'Starting buffer loop ', i+1
+    print('Starting buffer loop {}'.format(i+1))
     keepCellMaskNew = np.copy(keepCellMask)  # make a copy to edit that can be edited without changing the original
     ind = np.nonzero(keepCellMask == 0)[0]
     for i in range(len(ind)):
        iCell = ind[i]
        keepCellMaskNew[iCell] = keepCellMask[cellsOnCell[iCell,:nEdgesOnCell[iCell]]-1].max() # if any neighbor has a value of 1, then 1 will get assigned to iCell.
     keepCellMask = np.copy(keepCellMaskNew)  # after we've looped over all cells assign the new mask to the variable we need (either for another loop around the domain or to write out)
-    print '  Num of cells to keep:', sum(keepCellMask)
+    print('  Num of cells to keep: {}'.format(sum(keepCellMask)))
 
   # Now convert the keepCellMask to the cullMask
   cullCell[:] = np.absolute(keepCellMask[:]-1)  # Flip the mask for which ones to cull
@@ -100,13 +103,13 @@
 #  remove cells beyond a certain distance of ice extent
 elif maskmethod == 'distance':
 
-  print "Method: remove cells beyond a certain distance from existing ice"
+  print("Method: remove cells beyond a certain distance from existing ice")
 
   if thicknessMissing:
      sys.exit("Unable to perform 'numCells' method because thickness field was missing.")
 
   dist=float(options.distance)
-  print "Using a buffer distance of {} km".format(dist)
+  print("Using a buffer distance of {} km".format(dist))
   dist = dist * 1000.0 # convert to m
 
   keepCellMask = np.copy(cullCell[:])
@@ -118,7 +121,7 @@
 
   # mark the cells with ice first
   keepCellMask[thickness > 0.0] = 1
-  print 'Num of cells with ice:', sum(keepCellMask)
+  print('Num of cells with ice: {}'.format(sum(keepCellMask)))
 
   # find list of margin cells
   iceCells = np.nonzero(keepCellMask == 1)[0]
@@ -138,7 +141,7 @@
       ind = np.nonzero(((xCell-xCell[iCell])**2 + (yCell-yCell[iCell])**2)**0.5 < dist)[0]
       keepCellMask[ind] = 1
 
-  print '  Num of cells to keep:', sum(keepCellMask)
+  print('  Num of cells to keep:'.format(sum(keepCellMask)))
 
   # Now convert the keepCellMask to the cullMask
   cullCell[:] = np.absolute(keepCellMask[:]-1)  # Flip the mask for which ones to cull
@@ -147,14 +150,14 @@
 # =========
 #  cut out beyond some radius (good for the dome)
 elif maskmethod == 'radius':
-  print "Method: remove cells beyond a radius"
+  print("Method: remove cells beyond a radius")
   ind = np.nonzero( (xCell[:]**2 + yCell[:]**2)**0.5 > 26000.0 )
   cullCell[ind] = 1
 
 # =========
 #  cut off some fraction of the height/width on all 4 sides - useful for cleaning up a mesh from periodic_general
 elif maskmethod == 'edgeFraction':
-  print "Method: remove a fraction from all 4 edges"
+  print("Method: remove a fraction from all 4 edges")
   frac=0.025
 
   cullCell[:] = 0
@@ -175,7 +178,7 @@
 # =========
 
 
-print 'Num of cells to cull:', sum(cullCell[:])
+print('Num of cells to cull: {}'.format(sum(cullCell[:])))
 
 # =========
 # Try to add the new variable
@@ -205,6 +208,6 @@
   plt.show()
 
 f.close()
-print "cullMask generation complete."
+print("cullMask generation complete.")