<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Reading netCDF data
requires numpy and netCDF/HDF5 C libraries.
Github site
Step1:  Create a netCDF4.Dataset object
f is a Dataset object, representing an open netCDF file.
printing the object gives you summary information, similar to ncdump -h.
Step2:  Access a netCDF variable
variable objects stored by name in variables dict.
print the variable yields summary info (including all the attributes).
no actual data read yet (just have a reference to the variable object with metadata).
Step3:  List the Dimensions
All variables in a netCDF file have an associated shape, specified by a list of dimensions.
Let's list all the dimensions in this netCDF file.
Note that the MT dimension is special (unlimited), which means it can be appended to.
Step4:  Each variable has a dimensions and a shape attribute.
Step5:  Each dimension typically has a variable associated with it (called a coordinate variable).
Coordinate variables are 1D variables that have the same name as dimensions.
Coordinate variables and auxiliary coordinate variables (named by the coordinates attribute) locate values in time and space.
Step6:  Accessing data from a netCDF variable object
netCDF variables objects behave much like numpy arrays.
slicing a netCDF variable object returns a numpy array with the data.
Boolean array and integer sequence indexing behaves differently for netCDF variables than for numpy arrays. Only 1-d boolean arrays and integer sequences are allowed, and these indices work independently along each dimension (similar to the way vector subscripts work in fortran).
Step7:  What is the sea surface temperature and salinity at 50N, 140W?
Finding the latitude and longitude indices of 50N, 140W
The X and Y dimensions don't look like longitudes and latitudes
Use the auxilary coordinate variables named in the coordinates variable attribute, Latitude and Longitude
Step8:  Aha!  So we need to find array indices iy and ix such that Latitude[iy, ix] is close to 50.0 and Longitude[iy, ix] is close to -140.0 ...
Step9:  Now we have all the information we need to find our answer.
|----------+--------|
| Variable |  Index |
|----------+--------|
| MT       |      0 |
| Depth    |      0 |
| Y        | iy_min |
| X        | ix_min |
|----------+--------|
What is the sea surface temperature and salinity at the specified point?
Step10:  Remote data access via openDAP
Remote data can be accessed seamlessly with the netcdf4-python API
Access happens via the DAP protocol and DAP servers, such as TDS.
many formats supported, like GRIB, are supported "under the hood".
The following example showcases some nice netCDF features
Step11:  Missing values
when data == var.missing_value somewhere, a masked array is returned.
illustrate with soil moisture data (only defined over land)
white areas on plot are masked values over water.
Step12:  Packed integer data
There is a similar feature for variables with scale_factor and add_offset attributes.
short integer data will automatically be returned as float data, with the scale and offset applied.  
Dealing with dates and times
time variables usually measure relative to a fixed date using a certain calendar, with units specified like hours since YY
Step13:  Get index associated with a specified date, extract forecast data for that date.
Step14:  Get temp forecast for Boulder (near 40N, -105W)
use function getcloses_ij we created before...
Step15:  Simple multi-file aggregation
What if you have a bunch of netcdf files, each with data for a different year, and you want to access all the data as if it were in one file?
Step16:  MFDataset uses file globbing to patch together all the files into one big Dataset.
You can also pass it a list of specific files.
Limitations
Step17:  Closing your netCDF file
It's good to close netCDF files, but not actually necessary when Dataset is open for read access only.
<ASSISTANT_TASK:>
Python Code:
import netCDF4
import numpy as np
Explanation: Reading netCDF data
requires numpy and netCDF/HDF5 C libraries.
Github site: https://github.com/Unidata/netcdf4-python
Online docs: http://unidata.github.io/netcdf4-python/
Based on Konrad Hinsen's old Scientific.IO.NetCDF API, with lots of added netcdf version 4 features.
Developed by Jeff Whitaker at NOAA, with many contributions from users.
Interactively exploring a netCDF File
Let's explore a netCDF file from the Atlantic Real-Time Ocean Forecast System
first, import netcdf4-python and numpy
End of explanation
f = netCDF4.Dataset('data/rtofs_glo_3dz_f006_6hrly_reg3.nc')
print(f) 
Explanation: Create a netCDF4.Dataset object
f is a Dataset object, representing an open netCDF file.
printing the object gives you summary information, similar to ncdump -h.
End of explanation
print(f.variables.keys()) # get all variable names
temp = f.variables['temperature']  # temperature variable
print(temp) 
Explanation: Access a netCDF variable
variable objects stored by name in variables dict.
print the variable yields summary info (including all the attributes).
no actual data read yet (just have a reference to the variable object with metadata).
End of explanation
for d in f.dimensions.items():
    print(d)
Explanation: List the Dimensions
All variables in a netCDF file have an associated shape, specified by a list of dimensions.
Let's list all the dimensions in this netCDF file.
Note that the MT dimension is special (unlimited), which means it can be appended to.
End of explanation
temp.dimensions
temp.shape
Explanation: Each variable has a dimensions and a shape attribute.
End of explanation
mt = f.variables['MT']
depth = f.variables['Depth']
x,y = f.variables['X'], f.variables['Y']
print(mt)
print(x)                 
Explanation: Each dimension typically has a variable associated with it (called a coordinate variable).
Coordinate variables are 1D variables that have the same name as dimensions.
Coordinate variables and auxiliary coordinate variables (named by the coordinates attribute) locate values in time and space.
End of explanation
time = mt[:]  # Reads the netCDF variable MT, array of one element
print(time) 
dpth = depth[:] # examine depth array
print(dpth) 
xx,yy = x[:],y[:]
print('shape of temp variable: %s' % repr(temp.shape))
tempslice = temp[0, dpth > 400, yy > yy.max()/2, xx > xx.max()/2]
print('shape of temp slice: %s' % repr(tempslice.shape))
Explanation: Accessing data from a netCDF variable object
netCDF variables objects behave much like numpy arrays.
slicing a netCDF variable object returns a numpy array with the data.
Boolean array and integer sequence indexing behaves differently for netCDF variables than for numpy arrays. Only 1-d boolean arrays and integer sequences are allowed, and these indices work independently along each dimension (similar to the way vector subscripts work in fortran).
End of explanation
lat, lon = f.variables['Latitude'], f.variables['Longitude']
print(lat)
Explanation: What is the sea surface temperature and salinity at 50N, 140W?
Finding the latitude and longitude indices of 50N, 140W
The X and Y dimensions don't look like longitudes and latitudes
Use the auxilary coordinate variables named in the coordinates variable attribute, Latitude and Longitude
End of explanation
# extract lat/lon values (in degrees) to numpy arrays
latvals = lat[:]; lonvals = lon[:] 
# a function to find the index of the point closest pt
# (in squared distance) to give lat/lon value.
def getclosest_ij(lats,lons,latpt,lonpt):
    # find squared distance of every point on grid
    dist_sq = (lats-latpt)**2 + (lons-lonpt)**2  
    # 1D index of minimum dist_sq element
    minindex_flattened = dist_sq.argmin()    
    # Get 2D index for latvals and lonvals arrays from 1D index
    return np.unravel_index(minindex_flattened, lats.shape)
iy_min, ix_min = getclosest_ij(latvals, lonvals, 50., -140)
Explanation: Aha!  So we need to find array indices iy and ix such that Latitude[iy, ix] is close to 50.0 and Longitude[iy, ix] is close to -140.0 ...
End of explanation
sal = f.variables['salinity']
# Read values out of the netCDF file for temperature and salinity
print('%7.4f %s' % (temp[0,0,iy_min,ix_min], temp.units))
print('%7.4f %s' % (sal[0,0,iy_min,ix_min], sal.units))
Explanation: Now we have all the information we need to find our answer.
|----------+--------|
| Variable |  Index |
|----------+--------|
| MT       |      0 |
| Depth    |      0 |
| Y        | iy_min |
| X        | ix_min |
|----------+--------|
What is the sea surface temperature and salinity at the specified point?
End of explanation
import datetime
date = datetime.datetime.now()
# build URL for latest synoptic analysis time
URL = 'http://thredds.ucar.edu/thredds/dodsC/grib/NCEP/GFS/Global_0p5deg/GFS_Global_0p5deg_%04i%02i%02i_%02i%02i.grib2/GC' %\
(date.year,date.month,date.day,6*(date.hour//6),0)
# keep moving back 6 hours until a valid URL found
validURL = False; ncount = 0
while (not validURL and ncount < 10):
    print(URL)
    try:
        gfs = netCDF4.Dataset(URL)
        validURL = True
    except RuntimeError:
        date -= datetime.timedelta(hours=6)
        ncount += 1       
# Look at metadata for a specific variable
# gfs.variables.keys() will show all available variables.
sfctmp = gfs.variables['Temperature_surface']
# get info about sfctmp
print(sfctmp)
# print coord vars associated with this variable
for dname in sfctmp.dimensions:   
    print(gfs.variables[dname])
Explanation: Remote data access via openDAP
Remote data can be accessed seamlessly with the netcdf4-python API
Access happens via the DAP protocol and DAP servers, such as TDS.
many formats supported, like GRIB, are supported "under the hood".
The following example showcases some nice netCDF features:
We are seamlessly accessing remote data, from a TDS server.
We are seamlessly accessing GRIB2 data, as if it were netCDF data.
We are generating metadata on-the-fly.
End of explanation
soilmvar = gfs.variables['Volumetric_Soil_Moisture_Content_depth_below_surface_layer']
# flip the data in latitude so North Hemisphere is up on the plot
soilm = soilmvar[0,0,::-1,:] 
print('shape=%s, type=%s, missing_value=%s' % \
      (soilm.shape, type(soilm), soilmvar.missing_value))
import matplotlib.pyplot as plt
%matplotlib inline
cs = plt.contourf(soilm)
Explanation: Missing values
when data == var.missing_value somewhere, a masked array is returned.
illustrate with soil moisture data (only defined over land)
white areas on plot are masked values over water.
End of explanation
from netCDF4 import num2date, date2num, date2index
timedim = sfctmp.dimensions[0] # time dim name
print('name of time dimension = %s' % timedim)
times = gfs.variables[timedim] # time coord var
print('units = %s, values = %s' % (times.units, times[:]))
dates = num2date(times[:], times.units)
print([date.strftime('%Y-%m-%d %H:%M:%S') for date in dates[:10]]) # print only first ten...
Explanation: Packed integer data
There is a similar feature for variables with scale_factor and add_offset attributes.
short integer data will automatically be returned as float data, with the scale and offset applied.  
Dealing with dates and times
time variables usually measure relative to a fixed date using a certain calendar, with units specified like hours since YY:MM:DD hh-mm-ss.
num2date and date2num convenience functions provided to convert between these numeric time coordinates and handy python datetime instances.  
date2index finds the time index corresponding to a datetime instance.
End of explanation
from datetime import datetime, timedelta
date = datetime.now() + timedelta(days=3)
print(date)
ntime = date2index(date,times,select='nearest')
print('index = %s, date = %s' % (ntime, dates[ntime]))
Explanation: Get index associated with a specified date, extract forecast data for that date.
End of explanation
lats, lons = gfs.variables['lat'][:], gfs.variables['lon'][:]
# lats, lons are 1-d. Make them 2-d using numpy.meshgrid.
lons, lats = np.meshgrid(lons,lats)
j, i = getclosest_ij(lats,lons,40,-105)
fcst_temp = sfctmp[ntime,j,i]
print('Boulder forecast valid at %s UTC = %5.1f %s' % \
      (dates[ntime],fcst_temp,sfctmp.units))
Explanation: Get temp forecast for Boulder (near 40N, -105W)
use function getcloses_ij we created before...
End of explanation
!ls -l data/prmsl*nc
Explanation: Simple multi-file aggregation
What if you have a bunch of netcdf files, each with data for a different year, and you want to access all the data as if it were in one file?
End of explanation
mf = netCDF4.MFDataset('data/prmsl*nc')
times = mf.variables['time']
dates = num2date(times[:],times.units)
print('starting date = %s' % dates[0])
print('ending date = %s'% dates[-1])
prmsl = mf.variables['prmsl']
print('times shape = %s' % times.shape)
print('prmsl dimensions = %s, prmsl shape = %s' %\
     (prmsl.dimensions, prmsl.shape))
Explanation: MFDataset uses file globbing to patch together all the files into one big Dataset.
You can also pass it a list of specific files.
Limitations:
It can only  aggregate the data along the leftmost dimension of each variable.
only works with NETCDF3, or NETCDF4_CLASSIC formatted files.
kind of slow.
End of explanation
f.close()
gfs.close()
Explanation: Closing your netCDF file
It's good to close netCDF files, but not actually necessary when Dataset is open for read access only.
End of explanation
<END_TASK>