NGOM Interpolation
The NGOM step 2 model 6 - i.e. MODAS nowcast from Rich Patchen - files were received containing a orthogonal curvilinear grid rather than the standard regular grid requested.
A Python/Numpy script was written to interpolate all the fields in those files to the standard regular grid.
The interpolated files are approximately three times larger (6.5 Gb) than the original files (2.3 Gb), and the interpolation procedure requires about two hours for each file.
The script follows:
#!/usr/bin/python2.5
# encoding: utf-8
#
"""
gom.py
Created on 2011-09-25
Copyright (c) 2009 Steven K. Baum
Version 1:
"""
__docformat__ = "restructuredtext en"
import os, sys, fnmatch
import netCDF4
#import time,datetime
import matplotlib
from octant import csa
from netCDF4 import num2date, date2num
from mpl_toolkits.basemap import Basemap
import numpy as np
#from time import strftime
# Provide names for input and output files.
infile = netCDF4.Dataset('GOMHF-2011.3.101.GOMEX-PPP.nc','r+')
outfile = netCDF4.Dataset('test.nc','w',format='NETCDF3_CLASSIC')
#outfile = netCDF4.Dataset('GOMHF-2011.3.101.GOMEX-PPP-rg.nc','w',format='NETCDF3_CLASSIC')
csabathy = "/usr/local/bin/csabathy"
lonmin = -98.
lonmax = -78.
latmin = 18.
latmax = 31.0
# Set projection parameters.
lonzero = (lonmax - lonmin)/2.
latzero = (latmax - latmin)/2.
proj = Basemap(projection='lcc',resolution='i',llcrnrlon=lonmin,llcrnrlat= latmin,
urcrnrlon=lonmax,urcrnrlat=latmax,lat_0=latzero,lon_0=lonzero)
# lon - -98 to -78, dlon = 0.5
# lat - 18 to 31, dlat = .05
nlon = 401
nlat = 261
times = infile.variables['time'][:]
tstep = len(infile.variables['time'])
nsigma = len(infile.variables['sigma'])
print " nsigma = ",nsigma
sigmas = infile.variables['sigma'][:]
lon_2d_src = infile.variables['lon'][:]
lat_2d_src = infile.variables['lat'][:]
# Flatten the 2-D arrays.
lon_1d_src = lon_2d_src[:,:].ravel()
lat_1d_src = lat_2d_src[:,:].ravel()
# Convert lon/lat to x/y using proj
x_src,y_src = proj(lon_1d_src,lat_1d_src)
# Create lons/lats for regular array.
lon_dst = np.arange(-98,-77.95,.05)
lat_dst = np.arange(18,31.05,.05)
# Create 2-D arrays from the 1-D arrays using vstack.
lon_2d_dst = np.vstack((nlat*[lon_dst[:]]))
lat_2d_dst = np.vstack((nlon*[lat_dst[:]]))
lat_2d_dst = np.swapaxes(lat_2d_dst,0,1)
# Flatten the 2-D arrays for csabathy.
lon_1d_dst = lon_2d_dst[:,:].ravel()
lat_1d_dst = lat_2d_dst[:,:].ravel()
# Convert lon/lat to x/y using proj
x_dst,y_dst = proj(lon_1d_dst,lat_1d_dst)
outfile.comment = "Regridded file GOMHF-2011.3.101.GOMEX-PPP.nc to the GOMEX standard grid."
outfile.comment = "See http://stommel.tamu.edu/~baum/gomex-regrid.html for details."
# Create dimensions for output fields.
lon = outfile.createDimension('lon',nlon)
lat = outfile.createDimension('lat',nlat)
outfile.createDimension('sigma',nsigma)
outfile.createDimension('time',None)
# Create variables.
lon = outfile.createVariable('lon','f4',('lon',))
lon.long_name = "Longitude"
lon.units = "degrees_east"
lon.standard_name = "longitude"
lat = outfile.createVariable('lat','f4',('lat',))
lat.long_name = "Latitude"
lat.units = "degrees_north"
lat.standard_name = "latitude"
u = outfile.createVariable('u','f4',('time','sigma','lat','lon',))
u.long_name = "Eastward Water Velocity"
u.units = "m/s"
u.standard_name = "eastward_sea_water_velocity"
v = outfile.createVariable('v','f4',('time','sigma','lat','lon',))
v.long_name = "Northward Water Velocity"
v.units = "m/s"
v.standard_name = "northward_sea_water_velocity"
w = outfile.createVariable('w','f4',('time','sigma','lat','lon',))
w.long_name = "Vertical Water Velocity"
w.units = "m/s"
w.standard_name = "upward_sea_water_velocity"
temp = outfile.createVariable('temp','f4',('time','sigma','lat','lon',))
temp.long_name = "Temperature"
temp.units = "Celsius"
temp.standard_name = "sea_water_temperature"
salt = outfile.createVariable('salt','f4',('time','sigma','lat','lon',))
salt.long_name = "Salinity"
salt.units = "ppt"
salt.standard_name = "sea_water_salinity"
sigma = outfile.createVariable('sigma','f4',('sigma',))
sigma.long_name = "Sigma Stretched Vertical Coordinate at Nodes"
sigma.units = "sigma_level"
sigma.positive = "down"
sigma.standard_name = "ocean_sigma_coordinate"
sigma.formula_terms = "sigma: sigma eta: zeta depth: depth"
zeta = outfile.createVariable('zeta','f4',('time','lat','lon',))
zeta.long_name = "Water Surface Elevation"
zeta.units = "meters"
zeta.positive = "up"
zeta.standard_name = "sea_surface_elevation"
timeout = outfile.createVariable('time','f4',('time',))
timeout.long_name = "Time"
timeout.base_date = 2000, 8, 1, 0
timeout.units = "days since 2000-08-01 0:00:00 00:00"
timeout.standard_name = "time"
# Write the new lon/lat fields to the output file.
lon[:] = lon_dst
lat[:] = lat_dst
sigma[:] = sigmas
#tstep = 2
# Loop over number of time fields.
# ===========================================================================
for ts in range(0, tstep - 1, 1):
# ===========================================================================
print " ts = ",ts," times[ts] = ",times[ts]
timeout[ts] = times[ts]
# Extract 2-D zeta slice.
zeta_2d_src = infile.variables['zeta'][ts,:,:]
# print " zeta[3,3] = ",zeta_2d_src[50,50]
# Flatten 2-D zeta to 1-D zeta.
zeta_1d_src = zeta_2d_src[:,:].ravel()
# Interpolate the sucker.
csa_interp = csa.CSA(x_src,y_src,zeta_1d_src)
zeta_1d_dst = csa_interp(x_dst,y_dst)
# Reshape the U and V fields into 2D.
zeta_2d_dst = zeta_1d_dst.reshape(nlat,nlon)
# Write the interpolated zeta field for this time step to the output file.
zeta[ts,:] = zeta_2d_dst
# Loop over the freaking levels.
for sigma in range(0, nsigma, 1):
print " ts = ",ts," sigma = ",sigma
# Extract velocity fields.
print " Interpolating and writing u..."
u_2d_src = infile.variables['u'][ts,sigma,:,:]
# print " u_2d_src = ",u_2d_src[50,50]
u_1d_src = u_2d_src[:,:].ravel()
# print " u_1d_src = ",u_1d_src[2500]
csa_interp = csa.CSA(x_src,y_src,u_1d_src)
u_1d_dst = csa_interp(x_dst,y_dst)
u_2d_dst = u_1d_dst.reshape(nlat,nlon)
u[ts,sigma,:,:] = u_2d_dst
print " Interpolating and writing v..."
v_2d_src = infile.variables['v'][ts,sigma,:,:]
v_1d_src = v_2d_src[:,:].ravel()
csa_interp = csa.CSA(x_src,y_src,v_1d_src)
v_1d_dst = csa_interp(x_dst,y_dst)
v_2d_dst = v_1d_dst.reshape(nlat,nlon)
v[ts,sigma,:,:] = v_2d_dst
print " Interpolating and writing w..."
w_2d_src = infile.variables['w'][ts,sigma,:,:]
w_1d_src = w_2d_src[:,:].ravel()
csa_interp = csa.CSA(x_src,y_src,w_1d_src)
w_1d_dst = csa_interp(x_dst,y_dst)
w_2d_dst = w_1d_dst.reshape(nlat,nlon)
w[ts,sigma,:,:] = w_2d_dst
print " Interpolating and writing temp..."
temp_2d_src = infile.variables['temp'][ts,sigma,:,:]
temp_1d_src = temp_2d_src[:,:].ravel()
csa_interp = csa.CSA(x_src,y_src,temp_1d_src)
temp_1d_dst = csa_interp(x_dst,y_dst)
temp_2d_dst = temp_1d_dst.reshape(nlat,nlon)
temp[ts,sigma,:,:] = temp_2d_dst
print " Interpolating and writing salt..."
salt_2d_src = infile.variables['salt'][ts,sigma,:,:]
salt_1d_src = salt_2d_src[:,:].ravel()
csa_interp = csa.CSA(x_src,y_src,salt_1d_src)
salt_1d_dst = csa_interp(x_dst,y_dst)
salt_2d_dst = salt_1d_dst.reshape(nlat,nlon)
salt[ts,sigma,:,:] = salt_2d_dst
infile.close()
outfile.close()