<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 The netCDF file format
popular scientific file format for ocean and atmospere gridded datasets
netCDF is a collection of formats for storing arrays
netCDF classic
more widespread
2 GB file limit
often preffered for distributing products
netCDF 64 bit offset
supports larger files
NetCDF4
based on HDF5
compression
multiple unlimited variables
new types inc. user defined
herarchical groups
Developed by Unidata-UCAR with the aim of storing climate model data (3D+time)
Auxilary information about each variable can be added
Readable text equivalent called CDL (use ncdump/ncgen)
Can be used with Climate and Forecast (CF) data convention
http
Step1:  The main advantages of using xarray versus plain netCDF4 are
Step2:  ...or import local dataset
Step3:  Extract variable from dataset
Step4:  Access variable attributes
Step5:  Accessing data values
Step6:  Indexing and selecting data
From http
Step7:  Define selection using nearest value
Step8:  Plotting
Step9:  Arithmetic operations
Step10:  Calculate average along a dimension
Step11:  A dataset can easily be saved to a netCDF file
Step12:  Exercise
Extract the bathymetry
Extract the seabed temperature    isel(level=0)
Produce a scatter plot of depth vs. seabed temperature
<ASSISTANT_TASK:>
Python Code:
# Import everything that we are going to need... but not more
import pandas as pd
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap, cm
%matplotlib inline
DF=pd.DataFrame.from_items([('A', [1, 2, 3]), ('B', [4, 5, 6])],                         
                        orient='index', columns=['one', 'two', 'three'])
DF.mean(0)
pd_s=pd.Series(range(3), index=list('abc'), name='foo')
print(pd_s)
print()
#conver 1D series to ND aware dataArray 
print(xr.DataArray(pd_s))
Explanation: The netCDF file format
popular scientific file format for ocean and atmospere gridded datasets
netCDF is a collection of formats for storing arrays
netCDF classic
more widespread
2 GB file limit
often preffered for distributing products
netCDF 64 bit offset
supports larger files
NetCDF4
based on HDF5
compression
multiple unlimited variables
new types inc. user defined
herarchical groups
Developed by Unidata-UCAR with the aim of storing climate model data (3D+time)
Auxilary information about each variable can be added
Readable text equivalent called CDL (use ncdump/ncgen)
Can be used with Climate and Forecast (CF) data convention
http://cfconventions.org/
Data model:
Dimensions:describe the axes of the data arrays.
Variables: N-dimensional arrays of data.
Attributes: annotate variables or files with small notes or supplementary metadata.
Example for an ocean model dataset:
Dimensions
lat
lon
depth
time
Variable
Temperature
Salinity
Global Attibutes
Geographic grid type
History
Variable attributes (Temperature)
Long_name: "sea water temperature" 
Missing_value: 1.09009E36
Units: deg. C
range: -2:50
Tools for working with netCDF files
C and Fortran libraries
Used to underpin interfaces to other languages such as python (e.g. python package netCDF4)
Include ncdump/ncgen to convert to and from human readable CDL format.
nco tools http://nco.sourceforge.net/nco.html
Command line utilities to extract, create and operate data in netCDF files.
&gt; ncks -v u_wind -d lat,50.,51. -d lon,0.,5 inputfile.nc outputfile.nc
Viewers
ncview, pyncview, panoply, etc...
Imported by many packages
ArcGIS, QGIS, Surfer,Ferret etc...
Working With netCDF files using xarray
Alternative to plain netCDF4 access from python. 
Brings the power of pandas to environmental sciences, by providing N-dimensional variants of the core pandas data structures:
worth using for multidimensional data even when not using 
| Pandas | xarray  |
|---|---|
| Series  | DataArray  |
| DataFrame  | Dataset  |
DataArray uses names of dimensions making it easier to track than by using axis numbers. It is possible to write:
da.sel(time='2000-01-01') or da.mean(dim='time') 
intead of df.mean(0)
HTML documentation: http://xarray.pydata.org/
Example dataset
http://xray.readthedocs.io/en/stable/data-structures.html
<br>
<img src="../figures/dataset-diagram.png">
End of explanation
# Naughty datasets might require decode_cf=False
# Here it just needed decode_times=False
naughty_data = xr.open_dataset(
      'http://iridl.ldeo.columbia.edu/SOURCES/.OSU/.PRISM/.monthly/dods',
      decode_times=False)
naughty_data
Explanation: The main advantages of using xarray versus plain netCDF4 are:
intelligent selection along labelled dimensions (and also indexes)
groupby operations
data alignment
IO (netcdf)
conversion from and to Pandas.DataFrames
Lets peek inside our example file using ncdump
$ ncdump data/cefas_GETM_nwes.nc | more
netcdf cefas_GETM_nwes {
dimensions:
        latc = 360 ;
        lonc = 396 ;
        time = UNLIMITED ; // (6 currently)
        level = 5 ;
variables:
        double bathymetry(latc, lonc) ;
                bathymetry:units = "m" ;
                bathymetry:long_name = "bathymetry" ;
                bathymetry:valid_range = -5., 4000. ;
                bathymetry:_FillValue = -10. ;
                bathymetry:missing_value = -10. ;
        float h(time, level, latc, lonc) ;
                h:units = "m" ;
                h:long_name = "layer thickness" ;
                h:_FillValue = -9999.f ;
                h:missing_value = -9999.f ;
        double latc(latc) ;
                latc:units = "degrees_north" ;
        double level(level) ;
                level:units = "level" ;
        double lonc(lonc) ;
                lonc:units = "degrees_east" ;
        float temp(time, level, latc, lonc) ;
                temp:units = "degC" ;
                temp:long_name = "temperature" ;
                temp:valid_range = -2.f, 40.f ;
                temp:_FillValue = -9999.f ;
                temp:missing_value = -9999.f ;
        double time(time) ;
                time:long_name = "time" ;
                time:units = "seconds since 1996-01-01 00:00:00" ;
Now let's go back to python and use xarray.
Import remote dataset
xarray supports OpenDAP. This means that a dataset can be accessed remotely and subsetted as needed. Only the selected parts are downloaded.
End of explanation
GETM = xr.open_dataset('../data/cefas_GETM_nwes.nc4')
GETM
GETM.dims
print(type(GETM.coords['latc']))
GETM.coords['latc'].shape
# List name of dataset attributes
GETM.attrs.keys()
# List variable names
GETM.data_vars.keys()
Explanation: ...or import local dataset
End of explanation
temp=GETM['temp']
print(type( temp ))
temp.shape
Explanation: Extract variable from dataset
End of explanation
# print varaible attributes
for at in temp.attrs:
    print(at+':\t\t',end=' ')
    print(temp.attrs[at])
Explanation: Access variable attributes
End of explanation
temp[0,0,90,100]
Explanation: Accessing data values
End of explanation
#positional by integer
print( temp[0,2,:,:].shape )
# positional by label
print( temp.loc['1996-02-02T01:00:00',:,:,:].shape )
# by name and integer
print( temp.isel(level=1,latc=90,lonc=100).shape )
# by name and label
print( temp.sel(time='1996-02-02T01:00:00').shape )
#temp.loc
Explanation: Indexing and selecting data
From http://xarray.pydata.org/
<br>
<img src="../figures/xarray_indexing_table.png">
End of explanation
#GETM.sel(level=1)['temp']
GETM['temp'].sel(level=1,lonc=-5.0,latc=-50.0, method='nearest')
try:
    GETM['temp'].sel(level=1,lonc=-5.0,latc=-50.0, method='nearest',tolerance=0.5)
except KeyError:
    print('ERROR: outside tolerance of '+str(0.5))
Explanation: Define selection using nearest value
End of explanation
# Define a general mapping function using basemap
def do_map(var,title,units):
    latc=GETM.coords['latc'].values
    lonc=GETM.coords['lonc'].values
    # create figure and axes instances
    fig = plt.figure()
    ax = fig.add_axes()
    # create polar stereographic Basemap instance.
    m = Basemap(projection='stere', lon_0=0.,lat_0=60.,
            llcrnrlat=49,urcrnrlat=60,
            llcrnrlon=-10,urcrnrlon=15,resolution='l')
    # bondaries resolution can be 'c','l','i','h' or 'f'
    m.drawcoastlines(linewidth=0.5)
    m.fillcontinents(color='0.8')
    parallels = np.arange(-45,70,5)
    m.drawparallels(parallels,labels=[1,0,0,0],fontsize=10)
    m.drawparallels?
    meridians = np.arange(-15,20,5)
    m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10)
    # create arrays of coordinates for contourf
    lon2d,lat2d=np.meshgrid(lonc,latc)
    # draw filled contours.
    m.contourf(lon2d,lat2d,var,50,latlon=True)
    # add colorbar.
    cbar = m.colorbar(cmap=plt.cm.coolwarm,location='right')
    cbar.set_label(units)
    # add title
    plt.title(title)
# Extract attributes
units=GETM['temp'].attrs['units']
var_long_name=GETM['temp'].attrs['long_name']
# and plot
do_map(var=time_ave.sel(level=21),
       units=units,
       title='Time averaged '+var_long_name)
# But often, this will do
time_ave.sel(level=21).plot()
Explanation: Plotting
End of explanation
top=GETM['temp'].isel(time=0,level=4)
bottom=GETM['temp'].isel(time=0,level=0)
diff=top-bottom
diff.plot()
Explanation: Arithmetic operations
End of explanation
# average over time
time_ave = GETM['temp'].mean('time')
#average over time and level (vertical)
timelev_ave=GETM['temp'].mean(['time','level'])
timelev_ave.plot()
#zonal average (vertical)
timelon_ave=GETM['temp'].mean(['time','lonc']).isel(level=4)
timelon_ave.plot()
Explanation: Calculate average along a dimension
End of explanation
ds=GETM[['temp']].mean('time','level')
ds.to_netcdf('../data/temp_avg_level_time.nc')
print(type( GETM[['temp']]) )
print(type( GETM['temp'])   )
Explanation: A dataset can easily be saved to a netCDF file
End of explanation
# bathy = GETM
# bedtemp=GETM
# plt.scatter(  , ,marker='.')
Explanation: Exercise
Extract the bathymetry
Extract the seabed temperature    isel(level=0)
Produce a scatter plot of depth vs. seabed temperature
End of explanation
<END_TASK>