<p><b>mhoffman@lanl.gov</b> 2012-05-04 15:58:19 -0600 (Fri, 04 May 2012)</p><p>BRANCH COMMIT<br>
<br>
Created python script to copy the fields from a CISM netCDF file to a planar MPAS-format netCDF file (e.g. periodic_hex). Useful, e.g, for setting up a Greenland run. Also will provide a template for copying other regular grid files into MPAS format. Uses bilinear interpolation.<br>
</p><hr noshade><pre><font color="gray">Added: branches/land_ice_projects/grid_tools/copy_CISM_grid_to_MPAS_grid.py
===================================================================
--- branches/land_ice_projects/grid_tools/copy_CISM_grid_to_MPAS_grid.py (rev 0)
+++ branches/land_ice_projects/grid_tools/copy_CISM_grid_to_MPAS_grid.py 2012-05-04 21:58:19 UTC (rev 1865)
@@ -0,0 +1,206 @@
+#!/usr/bin/python
+# Copy fields from a regular CISM grid into a pre-existing MPAS grid
+
+import sys, numpy
+from Scientific.IO.NetCDF import *
+import math
+#import scipy.interpolate
+
+#----------------------------
+# Define which time levels you want to use in the two files! (0-based indexing in python)
+timelev = 0
+timelevout = 0
+#----------------------------
+
+
+
+# Check to see if the grid files were specified
+if len(sys.argv) == 3:
+ infilename = sys.argv[1]
+ outfilename = sys.argv[2]
+else:
+ print '</font>
<font color="blue">Usage: python *.py [CISM_GRID.NC] [MPAS_GRID.NC]'
+ print len(sys.argv)
+ sys.exit(0)
+
+#----------------------------
+# Define needed functions
+def BilinearInterp(x, y, Value, xCell, yCell):
+ # Calculate bilinear interpolation of Value field from x, y to new ValueCell field (return value) at xCell, yCell
+ # This assumes that x, y, Value are regular CISM style grids and xCell, yCell, ValueCell are 1-D unstructured MPAS style grids
+ try:
+ dx = x[1] - x[0]
+ dy = y[1] - y[0]
+ ValueCell = numpy.zeros(xCell.shape)
+ for i in range(len(xCell)):
+ # Calculate the CISM grid cell indices (these are the lower index)
+ xgrid = math.floor( xCell[i] / dx )
+ if xgrid >= len(x) - 1:
+ xgrid = len(x) - 2
+ ygrid = math.floor( yCell[i] / dy )
+ if ygrid >= len(y) - 1:
+ ygrid = len(y) - 2
+ #print xgrid, ygrid
+ ValueCell[i] = Value[ygrid,xgrid] * (x[xgrid+1] - xCell[i]) * (y[ygrid+1] - yCell[i]) / (dx * dy) + \
+ Value[ygrid+1,xgrid] * (x[xgrid+1] - xCell[i]) * (yCell[i] - y[ygrid]) / (dx * dy) + \
+ Value[ygrid,xgrid+1] * (xCell[i] - x[xgrid]) * (y[ygrid+1] - yCell[i]) / (dx * dy) + \
+ Value[ygrid+1,xgrid+1] * (xCell[i] - x[xgrid]) * (yCell[i] - y[ygrid]) / (dx * dy)
+ #print ValueCell.max(), ValueCell.min()
+ except:
+ 'error in BilinearInterp'
+ return ValueCell
+#----------------------------
+
+
+# Open the input file, get needed dimensions
+try:
+ infile = NetCDFFile(infilename,'r')
+
+ # Get the CISM dimensions if they exist
+ try:
+ level = infile.dimensions['level']
+ except:
+ print 'Input file is missing the dimension level. Might not be a problem.'
+
+# try:
+# x1 = infile.dimensions['x1']
+# y1 = infile.dimensions['y1']
+# except:
+# print 'Input file is missing the dimensions x1, y1. '
+
+# try:
+# x1 = infile.dimensions['x1']
+# y1 = infile.dimensions['y1']
+# except:
+# print 'Input file is missing the dimensions x0, y0. '
+
+ # Get CISM location variables if they exist
+ try:
+ x1 = infile.variables['x1'][:]
+ dx1 = x1[1] - x1[0]
+ print 'x1 min/max/dx:', x1.min(), x1.max(), dx1
+ y1 = infile.variables['y1'][:]
+ dy1 = y1[1] - y1[0]
+ print 'y1 min/max/dx:', y1.min(), y1.max(), dy1
+ except:
+ print 'Input file is missing x1 and/or y1. Might not be a problem.'
+
+ try:
+ x0 = infile.variables['x0'][:]
+ print 'x0 min/max:', x0.min(), x0.max()
+ y0 = infile.variables['y0'][:]
+ print 'y0 min/max:', y0.min(), y0.max()
+ except:
+ print 'Input file is missing x0 and/or y0. Might not be a problem.'
+
+except:
+ sys.exit('Error: The input file specified is either missing or lacking needed dimensions/variables.')
+
+
+
+# Open the output file, get needed dimensions & variables
+try:
+ outfile = NetCDFFile(outfilename,'r+')
+ nVertLevels = outfile.dimensions['nVertLevels']
+
+ # 1d vertical fields
+ layerThicknessFractions = outfile.variables['layerThicknessFractions'][:]
+
+ # '2d' spatial fields on cell centers
+ xCell = outfile.variables['xCell'][:]
+ print 'xCell min/max:', xCell.min(), xCell.max()
+ yCell = outfile.variables['yCell'][:]
+ print 'yCell min/max:', yCell.min(), yCell.max()
+
+except:
+ sys.exit('Error: The output grid file specified is either missing or lacking needed dimensions/variables.')
+
+
+#----------------------------
+# try each field. If it exists in the input file, it will be copied. If not, it will be skipped.
+# For now, include all variables defined for an MPAS land ice grid. If more are added (e.g. SMB) they will need to be added below.
+
+try:
+ thk = infile.variables['thk']
+ thickness = outfile.variables['thickness']
+ print '</font>
<font color="blue">thk min/max', thk[:].min(), thk[:].max()
+ thickness[timelevout,:] = BilinearInterp(x1, y1, thk[timelev,:,:], xCell, yCell)
+ print 'new thickness min/max', thickness[:].min(), thickness[:].max()
+ del thk, thickness
+except:
+ print '</font>
<font color="black">problem with thk field (e.g. not found in input file), skipping...</font>
<font color="blue">'
+
+
+try:
+ topg = infile.variables['topg']
+ bedTopography = outfile.variables['bedTopography']
+ print '</font>
<font color="blue">topg min/max', topg[:].min(), topg[:].max()
+ bedTopography[timelevout,:] = BilinearInterp(x1, y1, topg[timelev,:,:], xCell, yCell)
+ print 'new bedTopography min/max', bedTopography[:].min(), bedTopography[timelevout,:].max()
+ del topg, bedTopography
+except:
+ print '</font>
<font color="black">problem with topg field (e.g. not found in input file), skipping...</font>
<font color="blue">'
+
+
+try:
+ inbeta = infile.variables['beta']
+ beta = outfile.variables['beta']
+ print '</font>
<font color="blue">input beta min/max', inbeta[:].min(), inbeta[:].max()
+ beta[timelevout,:] = BilinearInterp(x0, y0, inbeta[timelev,:,:], xCell, yCell)
+ print 'new beta min/max', beta[:].min(), beta[:].max()
+ del inbeta, beta
+except:
+ print '</font>
<font color="black">problem with beta field (e.g. not found in input file), skipping...</font>
<font color="blue">'
+
+
+try:
+ temp = infile.variables['temp']
+ temperature = outfile.variables['temperature']
+ print '</font>
<font color="blue">input temp min/max', temp[:].min(), temp[:].max()
+ if level == nVertLevels + 2: # CISM includes the upper and lower boundaries as levels, MPAS does not.
+ #print range(1,level-1)
+ for i in range(1,level-1):
+ print 'Copying level ', i+1, ' of ', level , 'CISM levels (ignoring CISM b.c. temp levels)'
+ temperature[timelevout,:,i-1] = BilinearInterp(x1, y1, temp[timelev,i,:,:], xCell, yCell)
+ print 'new temperature min/max', temperature[:].min(), temperature[:].max()
+ else:
+ raise Exception
+ del temp, temperature
+except:
+ print '</font>
<font color="black">problem with temp field (e.g. not found in input file, differing number of layers), skipping...</font>
<font color="blue">'
+
+
+
+# Variables on edges...?
+#normalVelocity = outfile.variables['normalVelocity'][:]
+
+
+# Old implementation, including the super-slow scipy.interpolate object
+#try:
+# thk = infile.variables['thk'][:]
+# print 'thk min/max', thk.min(), thk.max()
+# #f = scipy.interpolate.RectBivariateSpline(y1, x1, thk[timelev, :, :], kx=1, ky=1)
+# #print 'got f'
+# #znew = f(yCell, xCell)
+# for i in range(len(xCell)):
+# # Calculate the CISM grid cell indices (these are the lower index)
+# xgrid = math.floor( xCell[i] / dx1 )
+# if xgrid >= len(x1)-1:
+# xgrid = len(x1) - 2
+# ygrid = math.floor( yCell[i] / dy1 )
+# if ygrid >= len(y1)-1:
+# ygrid = len(y1) - 2
+# #print xgrid, ygrid
+# thickness[0,i] = thk[timelev,ygrid,xgrid] * (x1[xgrid+1] - xCell[i]) * (y1[ygrid+1] - yCell[i]) / (dx1 * dy1) + \
+# thk[timelev,ygrid+1,xgrid] * (x1[xgrid+1] - xCell[i]) * (yCell[i] - y1[ygrid]) / (dx1 * dy1) + \
+# thk[timelev,ygrid,xgrid+1] * (xCell[i] - x1[xgrid]) * (y1[ygrid+1] - yCell[i]) / (dx1 * dy1) + \
+# thk[timelev,ygrid+1,xgrid+1] * (xCell[i] - x1[xgrid]) * (yCell[i] - y1[ygrid]) / (dx1 * dy1)
+# print 'new thickness min/max', thickness[:].min(), thickness[:].max()
+#except:
+# print 'problem with thk field (e.g. not found in input file), skipping...'
+
+
+infile.close()
+outfile.close()
+
+
</font>
</pre>