<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 In the first meeting of 2016 we had a brief look at the iris package API and one of its many capabilities. The data is in the same ERA-Interim reanalysis data.nc file, which was used in our previous exercises.
The idea is taken from python4oceanographers blog.
OK, Python, load the modules!
Step1:  Don't bother me with warnings!
Step2:  Read the NetCDF data file
The following switch is to automatically promote the reference surface variable associated with a dimensionless vertical coordinate. More info here
Step3:  Use the simplest loading method to open a NetCDF file as a iris.cube.Cubelist variable f.
Step4:  Let's take the horizontal wind components from this list
Step5:  Next, we construct a wind speed variable
Step6:  Note that we are performing arithmetic operations on an iris.cube.Cube, and the units of the derived variable are also correct. Another random example of the iris handling units conversions
Step7:  Subset the cube
iris makes it also easier to subset a cube in a human-readable way, without worrying about indices and instead using the physical limits. We will cut an area between 120 and 160 degrees longitude and 30 and 60 degrees latitude (Japan and surrounding area). It is done using iris.Constraint class and then extract() attribute of the original cube.
Step8:  iris trajectory
Import the trajectory submodule
Step9:  Define the starting and ending points of a cross-section
Step10:  Initialise a Trajectory with list of dictionaries pnts and set the number of sample positions to use along the trajectory to 100.
Step11:  Now the linearly spaced points of the trajectory are stored in its attribute sampled_points, which is a list of dictionaries. We could have got the same result just by applying numpy.linspace function, but this way is more convenient. However, we have to rearrange this data structure so that it can be fed to the interpolate() function
Step12:  We then can easily dissect our cube of wind speed along the defined trajectory.
Step13:  Note that the result is also an iris cube, whose dimensional coordinates that are not used in the interpolation remained preserved, as well as all other metadata. Also, the longitude and latitude coordinates moved to the auxiliary coordinate group.
Plot the result
Firstly, we need to extract the lon/lat arrays from the wspd_subset variable.
Step14:  We can also use the coordinates of the section to mark its location on a map.
Step15:  And then, we create a figure with two subplots
<ASSISTANT_TASK:>
Python Code:
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import iris
import matplotlib.pyplot as plt
import numpy as np
import os
%matplotlib inline
Explanation: In the first meeting of 2016 we had a brief look at the iris package API and one of its many capabilities. The data is in the same ERA-Interim reanalysis data.nc file, which was used in our previous exercises.
The idea is taken from python4oceanographers blog.
OK, Python, load the modules!
End of explanation
import warnings
warnings.filterwarnings('ignore')
Explanation: Don't bother me with warnings!
End of explanation
iris.FUTURE.netcdf_promote = True
Explanation: Read the NetCDF data file
The following switch is to automatically promote the reference surface variable associated with a dimensionless vertical coordinate. More info here: PR1159.
End of explanation
filepath = os.path.join(os.path.pardir,'data','data.nc') # ../data/data.nc on Unix
f = iris.load(filepath)
print(f)
Explanation: Use the simplest loading method to open a NetCDF file as a iris.cube.Cubelist variable f.
End of explanation
u, v = f[2], f[4]
print(u)
print('\nType of u variable: {}'.format(type(u)))
Explanation: Let's take the horizontal wind components from this list:
End of explanation
wspd = (u**2 + v**2)**0.5
wspd.rename('wind_magnitude')
print(wspd)
Explanation: Next, we construct a wind speed variable:
End of explanation
(123*u*u.coord('time')/u.coord('pressure_level')).units
Explanation: Note that we are performing arithmetic operations on an iris.cube.Cube, and the units of the derived variable are also correct. Another random example of the iris handling units conversions:
End of explanation
sublon = iris.Constraint(longitude=lambda cell: 120 <= cell <= 160)
sublat = iris.Constraint(latitude=lambda cell: 30 <= cell <= 60)
wspd_subset = wspd.extract(sublon & sublat)
print(wspd_subset)
import shapely.geometry as sgeom
fig, ax = plt.subplots(subplot_kw=dict(projection=ccrs.PlateCarree()))
ax.coastlines()
extent_box = sgeom.box(120,30,160,60)
ax.add_geometries([extent_box], ccrs.PlateCarree(), color='red', alpha=0.5, edgecolor='red', linewidth=2)
print('Subset location')
Explanation: Subset the cube
iris makes it also easier to subset a cube in a human-readable way, without worrying about indices and instead using the physical limits. We will cut an area between 120 and 160 degrees longitude and 30 and 60 degrees latitude (Japan and surrounding area). It is done using iris.Constraint class and then extract() attribute of the original cube.
End of explanation
from iris.analysis import trajectory
Explanation: iris trajectory
Import the trajectory submodule:
End of explanation
pnts = [{'longitude': 155, 'latitude': 35}, # start
        {'longitude': 125, 'latitude': 55} # end
       ]
Explanation: Define the starting and ending points of a cross-section:
End of explanation
traj = trajectory.Trajectory(pnts, sample_count=100)
traj
Explanation: Initialise a Trajectory with list of dictionaries pnts and set the number of sample positions to use along the trajectory to 100.
End of explanation
lon = [d['longitude'] for d in traj.sampled_points]
lat = [d['latitude'] for d in traj.sampled_points]
sampled_points = [('longitude', lon),
                  ('latitude', lat)]
Explanation: Now the linearly spaced points of the trajectory are stored in its attribute sampled_points, which is a list of dictionaries. We could have got the same result just by applying numpy.linspace function, but this way is more convenient. However, we have to rearrange this data structure so that it can be fed to the interpolate() function:
End of explanation
section = trajectory.interpolate(wspd_subset, sampled_points)
print(section)
Explanation: We then can easily dissect our cube of wind speed along the defined trajectory.
End of explanation
lon, lat = wspd_subset.coord('longitude').points, wspd_subset.coord('latitude').points
Explanation: Note that the result is also an iris cube, whose dimensional coordinates that are not used in the interpolation remained preserved, as well as all other metadata. Also, the longitude and latitude coordinates moved to the auxiliary coordinate group.
Plot the result
Firstly, we need to extract the lon/lat arrays from the wspd_subset variable.
End of explanation
seclon, seclat = section.coord('longitude').points, section.coord('latitude').points
Explanation: We can also use the coordinates of the section to mark its location on a map.
End of explanation
ilev = 0
it = 0
fig = plt.figure(figsize=(30,10))
ax1 = fig.add_subplot(121, projection=ccrs.PlateCarree())
ax1.coastlines('50m')
ax1.contourf(lon, lat, wspd_subset.data[it,ilev,...], cmap=plt.cm.viridis)
#--------Fancy formatting---------
gl = ax1.gridlines(crs=ccrs.PlateCarree(), # using the same projection
                  draw_labels=True, # add labels
                  linewidth=2, color='gray', alpha=0.5, linestyle='--') # grid line specs
# Remove labels above and on the right of the map (note that Python allows the double equality)
gl.xlabels_top = gl.ylabels_right = False
# Format the labels using the formatters imported from cartopy
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
#---------------------------------
ax1.plot(seclon, seclat, color='r', linestyle='', marker='o', linewidth=3)
ax1.set_title('Pressure level: {}hPa'.format(wspd_subset.coord('pressure_level').points[ilev]))
ax2 = fig.add_subplot(122)
ax2.contourf(np.arange(section.data.shape[2]), section.coord('pressure_level').points, section.data[it,...], cmap=plt.cm.viridis)
ax2.invert_yaxis()
ax2.set_ylabel('Pressure, hPa')
ax2.set_xlabel('Distance along the section, points')
fig.suptitle('Wind speed at {}'.format(wspd_subset.coord('time').units.num2date(wspd_subset.coord('time').points[0]).strftime('%Y-%b-%d')))
HTML(html)
Explanation: And then, we create a figure with two subplots: one with geoaxes for the map of the subset cube and one with pressure-distance axes for the cross-section.
Choose the first pressure level and the first time step:
End of explanation
<END_TASK>