The NOAA CO-OPS forecast models SJROFS uses the curvilinear grid model CH3D. The model grid was developed with a grid generation tool that allows regions of the grid to be undefined geographically, and therefore the lon/lat coordinate arrays contain missing values. CF Conventions (state "Missing values are not allowed in coordinate variables" do not allow coordinate arrays to contain missing values, and not surprisingly, the standard visualization tools that work with most curvilinear grids (e.g. godiva2) do not work with SJROFS grids.

But what if we could fill in the lon/lat values via interpolation so that we had no missing values? We explore this here, trying inverse distance interpolation. This unfortunately leads to non-monotonic behavior in lon/lat which also violates CF Conventions and leads to plotting problems. CF states that coordinate variables are a "numeric data type with values that are ordered monotonically".



In [1]:

    
from pylab import *
import netCDF4

rcParams['figure.figsize'] = 10,6



In [2]:

    
url='http://opendap.co-ops.nos.noaa.gov/thredds/dodsC/SJROFS/fmrc/Aggregated_7_day_SJROFS_Fields_Forecast_best.ncd'
nc = netCDF4.Dataset(url)
ncv = nc.variables



In [3]:

    
lon = ncv['lon'][:,:]
lat = ncv['lat'][:,:]
mask = ncv['mask'][:,:]
h = ncv['depth'][:,:]



In [4]:

    
print ncv['mask']









    



<type 'netCDF4.Variable'>
float32 mask(ny, nx)
    units: nondimensional
    long_name: Land Mask
    coordinates: lat lon 
unlimited dimensions: 
current shape = (105, 188)



In [5]:

    
print unique(mask)









    



[ 0.  5.]



In [22]:

    
lon3=ma.masked_where(mask<1,lon)
lat3=ma.masked_where(mask<1,lat)
plot(lon3,lat3,'k-',lon3.T,lat3.T,'k-');



In [7]:

    
plt.scatter(lon3, lat3, c=h, s=30, marker='o',edgecolors='none', cmap=mpl.cm.jet)









    Out[7]:





<matplotlib.collections.PathCollection at 0x46c6d90>

What do lon and lat look like in (i,j) space?



In [8]:

    
imshow(ma.masked_where(mask<1,lat))









    Out[8]:





<matplotlib.image.AxesImage at 0x35b3c10>



In [9]:

    
imshow(ma.masked_where(mask<1,lon))









    Out[9]:





<matplotlib.image.AxesImage at 0x35d9650>

try filling in the missing lon,lat values using inverse distance technique



In [10]:

    
(ny,nx) = shape(lon)
(ii,jj) = meshgrid(range(nx),range(ny))



In [11]:

    
iland = where(mask.ravel()==0)[0]
iwater = where(mask.ravel()==5)[0]



In [12]:

    
def distance_matrix(x0, y0, x1, y1):
    obs = np.vstack((x0, y0)).T
    interp = np.vstack((x1, y1)).T

    # Make a distance matrix between pairwise observations
    # Note: from <http://stackoverflow.com/questions/1871536>
    # (Yay for ufuncs!)
    d0 = np.subtract.outer(obs[:,0], interp[:,0])
    d1 = np.subtract.outer(obs[:,1], interp[:,1])

    return np.hypot(d0, d1)



In [13]:

    
def simple_idw(x, y, z, xi, yi):
    dist = distance_matrix(x,y, xi,yi)

    # In IDW, weights are 1 / distance
    weights = 1.0 / (dist*dist)

    # Make weights sum to one
    weights /= weights.sum(axis=0)

    # Multiply the weights for each interpolated point by all observed Z-values
    zi = np.dot(weights.T, z)
    return zi



In [14]:

    
lon2=lon.ravel()
lat2=lat.ravel()
jj = jj.ravel()
ii = ii.ravel()

Do the inverse distance interpolation in (i,j) space:



In [15]:

    
lon2[iland] = simple_idw(ii[iwater],jj[iwater],lon2[iwater],ii[iland],jj[iland])



In [16]:

    
lat2[iland] = simple_idw(ii[iwater],jj[iwater],lat2[iwater],ii[iland],jj[iland])

Well, how did our interpolation turn out?



In [17]:

    
lat2=lat2.reshape(shape(lat))
imshow(lat2)









    Out[17]:





<matplotlib.image.AxesImage at 0x493dfd0>



In [18]:

    
lon2=lon2.reshape(shape(lon))
imshow(lon2)









    Out[18]:





<matplotlib.image.AxesImage at 0x44ba410>

So what does the resulting mesh look like? Yeesh....



In [19]:

    
plot(lon2,lat2,'k-',lon2.T,lat2.T,'k-');



In [20]:

    
pcolormesh(lon2,lat2,h)









    Out[20]:





<matplotlib.collections.QuadMesh at 0x53f6210>



In [20]: