Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Land ice mesh tools to python 3 #249

Merged
merged 2 commits into from
May 8, 2019
Merged

Land ice mesh tools to python 3 #249

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
43895b1
Update print statement to function calls for python 3
xylar Feb 26, 2019
f00466a
Remove a tab (indentation should be with spaces)
xylar Feb 26, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
65 changes: 34 additions & 31 deletions landice/mesh_tools_li/create_landice_grid_from_generic_MPAS_grid.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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")
Expand All @@ -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')
Expand All @@ -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


# ============================================
Expand All @@ -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
Expand All @@ -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])
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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)

Expand All @@ -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'))
Expand All @@ -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'))
Expand All @@ -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'))
Expand All @@ -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'))
Expand All @@ -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[:])
Expand All @@ -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))
39 changes: 21 additions & 18 deletions landice/mesh_tools_li/define_cullMask.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand All @@ -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:
Expand All @@ -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.")

Expand All @@ -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
Expand All @@ -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
Expand All @@ -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[:])
Expand All @@ -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]
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -205,6 +208,6 @@
plt.show()

f.close()
print "cullMask generation complete."
print("cullMask generation complete.")


Loading