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)
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:
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
There are many ways of getting a gridded field from sparsely distributed observations. We will show two simple applications.
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)
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()
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()
In [12]:
print temperature_checked.min()
print temperature_checked.max()
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 [ ]: