Now that we know how to read the data, we will add some processing to generate a gridded field.
In the first part, we will use the loop on the files to store all the data, then in the second part we will compute a mean value for each grid point of the grid.

First we select a year and a month:


In [1]:
year = 2015
month = 7

In [2]:
%matplotlib inline
import glob
import os
import netCDF4
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from matplotlib  import colors

The directory where we store the data files:


In [3]:
basedir = "~/DataOceano/MyOcean/INSITU_GLO_NRT_OBSERVATIONS_013_030/monthly/" + str(year) + str(month).zfill(2) + '/'
basedir = os.path.expanduser(basedir)

Reading all the data

The loop is the same as the 2 previous examples, except that now we will keep the data in arrays that we create empty:


In [4]:
lon_checked, lat_checked, temperature_checked = np.array([]), np.array([]), np.array([])

filelist = sorted(glob.glob(basedir+'*.nc'))
k = 0
for datafiles in filelist:

    # print datafiles
    with netCDF4.Dataset(datafiles) as nc:
        
        lon = nc.variables['LONGITUDE'][:]
        lat = nc.variables['LATITUDE'][:]
        depth = nc.variables['DEPH'][:]
        POSITION_QC = nc.variables['POSITION_QC'][:]
        
        if depth.shape[1] == 1:

            try:
                TEMP_QC = nc.variables['TEMP_QC'][:, 0]
                temperature = nc.variables['TEMP'][:]         
                gooddata = np.where(np.logical_and((TEMP_QC == 1), (POSITION_QC == 1)))
               
                temperature = temperature[gooddata]
                
                temperature_checked = np.append(temperature_checked, temperature)
                lon_checked = np.append(lon_checked, lon[gooddata])
                lat_checked = np.append(lat_checked, lat[gooddata])
        
            except KeyError:
                k += 1
                # print 'No variable temperature in this file'

Now we have the coordinates and the temperature for all the files, we can save them in a file that we can re-use later.
To do so, we will use the numpy savetxt function. The function requires:

  • the name of the file where the data will be written,
  • the data to be saved, in the form of an array. To create a unique array from the different arrays (coordinates and temperature), we use the function c_.

In [5]:
datafile = './lon_lat_temperature_' + str(year) + '_' + str(month).zfill(2) + '.txt'
np.savetxt(datafile, np.c_[lon_checked, lat_checked, temperature_checked], fmt='%3.5f %3.5f %2.2f')

print 'Data saved in file ' + datafile


Data saved in file ./lon_lat_temperature_2015_07.txt

Creation of a gridded field

There are many ways of getting a gridded field from sparsely distributed observations. We will show two simple applications.

Linear interpolation

With scipy, the module interpolate provide many functions to perform interpolations.
In particular, griddata aims to interpolate unstructured D-dimensional data.


In [6]:
from scipy.interpolate import griddata

We need to specify the grid on which the observations have to be interpolated.
We construct a 5º by 5º grid between 70ºS and 70ºN latitude.


In [7]:
lon_interp, lat_interp = np.meshgrid(np.arange(-180, 180.5, 5.), np.arange(-70, 70., 5.))

The interpolated field is obtained as:


In [8]:
temperature_interp = griddata((lon_checked, lat_checked), temperature_checked, (lon_interp, lat_interp), method='linear')

After the interpolation and prior to the plot, it is necessary to mask the NaN values that could have been generated.


In [9]:
temperature_interp = np.ma.masked_where(np.isnan(temperature_interp), temperature_interp)

Plotting

Same code as in the previous examples. pcolormesh (pseudocolor plot of a 2-D array) is used for the representation of the gridded field.


In [10]:
fig = plt.figure(figsize=(12, 8))
tempmin, tempmax = 0., 30.
cmaptemp = plt.cm.RdYlBu_r
normtemp = colors.Normalize(vmin=tempmin, vmax=tempmax)
tempticks = np.arange(tempmin, tempmax+0.1,2.5)

m = Basemap(projection='moll', lon_0=0, resolution='c')
lon_interp_map, lat_interp_map = m(lon_interp, lat_interp)

pcm = m.pcolormesh(lon_interp_map, lat_interp_map, temperature_interp, cmap=cmaptemp, norm=normtemp)

cbar = plt.colorbar(pcm, extend='both', shrink=0.7)
cbar.set_label('$^{\circ}$C', rotation=0, ha='left')
m.drawcoastlines(linewidth=0.2)
m.fillcontinents(color = 'gray')
plt.title('Gridded temperature from surface drifters\n' + str(year) + '-' + str(month).zfill(2))
plt.show()


/usr/local/lib/python2.7/dist-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if self._edgecolors == str('face'):

The chosen method (linear) is far from the best, but the goal is simply to illustrate the tool with the drifter data.
We check that the min/max values are consistent with the data:


In [11]:
print temperature_interp.min()
print temperature_interp.max()


-2.95728333526
32.4052053629

In [12]:
print temperature_checked.min()
print temperature_checked.max()


-3.0
35.2999992371

Interpolation in a specific region

If we want to focus on a given area and increase the resolution, we just have to change the interpolation grid:


In [22]:
lonmin, lonmax, latmin, latmax = -25., 35., 20.0, 60.
deltalon, deltalat = 2., 2.

In [23]:
lon_interp, lat_interp = np.meshgrid(np.arange(lonmin, lonmax, deltalon), np.arange(latmin, latmax, deltalat))
temperature_interp = griddata((lon_checked, lat_checked), temperature_checked, (lon_interp, lat_interp), method='linear')
temperature_interp = np.ma.masked_where(np.isnan(temperature_interp), temperature_interp)

On the map we will also add the data locations.


In [24]:
fig = plt.figure(figsize=(12, 8))
m = Basemap(llcrnrlon=lonmin, llcrnrlat=latmin,
            urcrnrlon=lonmax, urcrnrlat=latmax, resolution='i')

lon_interp_map, lat_interp_map = m(lon_interp, lat_interp)
lon_checked_map, lat_checked_map = m(lon_checked, lat_checked)
pcm = m.pcolormesh(lon_interp_map, lat_interp_map, temperature_interp, cmap=cmaptemp, norm=normtemp)
plt.plot(lon_checked_map, lat_checked, 'ko', ms=0.1)

cbar = plt.colorbar(pcm, extend='both', shrink=0.9)
cbar.set_label('$^{\circ}$C', rotation=0, ha='left')
m.drawcoastlines(linewidth=0.2)
m.fillcontinents(color = 'gray')
plt.title('Gridded temperature from surface drifters\n' + str(year) + '-' + str(month).zfill(2), fontsize=20)
plt.show()


The figure illustrates the low number of available drifters in the Mediterranean Sea if compared to the Atlantic Ocean.


In [ ]: