Using python, openDAP, and matplotlib to plot CoRTAD SST on maps

SST values for the 4 tiles covering my model domain. The model grid is shown as red dots. Only every 10th grid cell is plotted.

In an earlier post I wrote on how to use e OPeNDAP to extract sea surface temperature (SST) data for a given region of the ocean. I have earlier used the CoRTAD version 3 data, but now there is a new and improved version 4 available which I wanted to use for my assimilation project in ROMS.  My plan is to use CoRTAD data to assimilate SST in the regional ocean modeling system (ROMS) using the IS4DVAR method. This requires me provide an observation file in NetCDF4 format that will be read by ROMS. Basically, the observation file contains all of the observations I want to assimilate into my ROMS model in sequential, time-dependent order. To create this file I therefore needed to write a script that can extract the SST values from the CoRTAD OPeNDAP server for my domain, interpolate the values to my ROMS grid, and finally write the values to file. Each observation needs to be accompanied by its longitude/latitude position, time of observation, observational value, and the corresponding position in the ROMS grid file. To start I define the base information for the CoRTAD OpenDAP server where all information is stored:

base="http://data.nodc.noaa.gov/thredds/dodsC/cortad/Version4/"

The area that I want to assimilate SST data is the North Sea. As the CoRTAD data is provided as tiles, I needed to combine 4 tiles to completely cover my model region. I extract the data for a given time step for four different tiles and combine them using this method:

def extractCORTADSST(name,t):
    """Routine that extracts the SST values for
    the specific tiles and time-period (t)"""
    cdf1,cdf2,cdf3,cdf4=openCoRTAD()

    filledSST1=cdf1.variables["FilledSST"][t,:,:]
    filledSST2=cdf2.variables["FilledSST"][t,:,:]
    filledSST3=cdf3.variables["FilledSST"][t,:,:]
    filledSST4=cdf4.variables["FilledSST"][t,:,:]

    offset=cdf1.variables["FilledSST"].__getattribute__('add_offset')
    cdf1.close();cdf2.close();cdf3.close();cdf4.close();
    filledMaskedSST1=filledSST1 - offset
    filledMaskedSST2=filledSST2 - offset
    filledMaskedSST3=filledSST3 - offset
    filledMaskedSST4=filledSST4 - offset

    """Now we have all the data in 4 different arrays that we need to concentate.
    First we add the horisontal tiles, and finally we stack the two horisontal ones on top
    of each other."""
    filledMaskedSST_lower=concatenate((filledMaskedSST1,filledMaskedSST2),axis=1)
    filledMaskedSST_upper=concatenate((filledMaskedSST3,filledMaskedSST4),axis=1)
    filledMaskedSST_all=concatenate((filledMaskedSST_upper,filledMaskedSST_lower),axis=0)

    """Flip the SST array to be consistent with order of latitude array"""
    filledMaskedSST_all=np.flipud(filledMaskedSST_all)

    """ Scale and offset is autmoatically detected and edited by netcdf, but
    we need to mask the values that are not filled."""
    filledMaskedSST_final=ma.masked_less(filledMaskedSST_all,-2.)

    print "Min and max of SST: %s - %s"%(filledMaskedSST_final.min(),filledMaskedSST_final.max())
    print "------------------------------\n"

    return filledMaskedSST_final

The arrays for the four tiles are combined into one array by the numpy function ‘concatenate’ Since the values in CoRTAD are in Kelvin we need to subtract the offset (273.15) to get the values in celsius. The values are actually also scaled by a factor of 0.01, but Python NetCDF4 interface takes care of that automatically for us. The time array in CoRTAD is defined as number of days since 1980.12.31 12:00:00 and the following method enables us to figure out what date a time index is equivalent to

def getCORTADtime():

    base="http://data.nodc.noaa.gov/thredds/dodsC/cortad/Version4/"
    file1="cortadv4_row01_col07.nc"
    filename1=base+file1
    cdf1=Dataset(filename1)

    print "Time: Extrating timedata from openDAP: %s"%(filename1)

    time=np.squeeze(cdf1.variables["time"][:])
    cdf1.close()
    return time

Call the time method on time index ‘t’ using:

time=getCortad.getCORTADtime()
refDate=datetime.datetime(1980,12,31,12,0,0)
currentDate=refDate + datetime.timedelta(days=float(time[t]))

After I had extracted the SST data for all tiles and stored the values into a large array, I interpolated the values to my ROMS grid. This grid is defined by my 2D longitude, latitude grid (grid_lon, grid_lat). To visualize my model domain relative to the full four tile CoRTAD domain, I extracted the boundary grid points from my grid and created a polygon.

def getPolygon(grid_lon,grid_lat):

    lon_bd = np.concatenate((grid_lon[:,0],grid_lon[-1,:],grid_lon[::-1,-1], grid_lon[0,::-1] ))
    lat_bd = np.concatenate((grid_lat[:,0],grid_lat[-1,:],grid_lat[::-1,-1], grid_lat[0,::-1] ))

    """ Save the polygon as array to plot later"""
    polygon_data=np.empty((2,len(lon_bd)))
    for k in xrange(len(lon_bd)):
        polygon_data[0,k]=lon_bd[k]
        polygon_data[1,k]=lat_bd[k]

    return polygon_data

This polygon array was then plotted onto my map using matplotlib

def drawGridBoundaries(ax,map,polygon_data):
    print "Plotting grid boundaries"
    polygon_data_xy=[]
    for i in xrange(len(polygon_data[0,:])):
        myx,myy=map(polygon_data[0,i],polygon_data[1,i])
        vertices=[myx,myy]
        polygon_data_xy.append(vertices)

    polygon_data_xy=np.asarray(polygon_data_xy)
    patch = Polygon(array(polygon_data_xy), facecolor='none',
        edgecolor=(.9, .3, .3, 1), linewidth=4)
    ax.add_patch(patch)

Finally, I plotted both the interpolated SST values and the full CoRTAD region using the regular matplotlib contourf function.I made the CoRTAD layer slightly translucent (alpha=0.5) to distinguish it from the interpolated SST values. The polygon boundary of my grid domain I plotted as a red line.

All off the python files can be downloaded here (zip file): CoRTAD (755)

Leave a Reply

Your email address will not be published. Required fields are marked *