<jupyter_start><jupyter_text>Continuing the previous week's topic, today we went through the basics of loading netCDF datasets. In this notebook we explore a netCDF file from the Atlantic Real-Time Ocean Forecast System (downloaded from the Unidata Python Workshop's [GitHub repository](https://unidata.github.io/unidata-python-workshop/). The filename is [rtofs_glo_3dz_f006_6hrly_reg3.nc](https://github.com/Unidata/tds-python-workshop/blob/master/data/rtofs_glo_3dz_f006_6hrly_reg3.nc).<jupyter_code>import os
path_to_file = os.path.join(os.path.pardir, 'data', 'rtofs_glo_3dz_f006_6hrly_reg3.nc')
print(path_to_file)<jupyter_output>../data/rtofs_glo_3dz_f006_6hrly_reg3.nc<jupyter_text>First, let's load essential modules<jupyter_code>from __future__ import division, print_function
import netCDF4 as nc
import numpy as np<jupyter_output><empty_output><jupyter_text>Create a netCDF4.Dataset object from file<jupyter_code>f = nc.Dataset(path_to_file)
print(f)<jupyter_output><class 'netCDF4._netCDF4.Dataset'>
root group (NETCDF4_CLASSIC data model, file format HDF5):
    Conventions: CF-1.0
    title: HYCOM ATLb2.00
    institution: National Centers for Environmental Prediction
    source: HYCOM archive file
    experiment: 90.9
    history: archv2ncdf3z
    dimensions(sizes): MT(1), Y(850), X(712), Depth(10)
    variables(dimensions): float64 [4mMT[0m(MT), float64 [4mDate[0m(MT), float32 [4mDepth[0m(Depth), int32 [4mY[0m(Y), int32 [4mX[0m(X), float32 [4mLatitude[0m(Y,X), float32 [4mLongitude[0m(Y,X), float32 [4mu[0m(MT,Depth,Y,X), float32 [4mv[0m(MT,Depth,Y,X), float32 [4mtemperature[0m(MT,Depth,Y,X), float32 [4msalinity[0m(MT,Depth,Y,X)
    groups:<jupyter_text>So what have we got here?- **`f`** is a `Dataset` object, representing an open netCDF file.- printing the object gives you summary information, similar to *`ncdump -h`*. You can also get a list of global attributes:<jupyter_code>f.ncattrs()<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>for d in f.dimensions.items():
    print(d)<jupyter_output>('MT', <class 'netCDF4._netCDF4.Dimension'> (unlimited): name = 'MT', size = 1
)
('Y', <class 'netCDF4._netCDF4.Dimension'>: name = 'Y', size = 850
)
('X', <class 'netCDF4._netCDF4.Dimension'>: name = 'X', size = 712
)
('Depth', <class 'netCDF4._netCDF4.Dimension'>: name = 'Depth', size = 10
)<jupyter_text>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).<jupyter_code>type(f.variables)
print(f.variables.keys()) # get all variable names

# or like this: [*f.variables.keys()]<jupyter_output>odict_keys(['MT', 'Date', 'Depth', 'Y', 'X', 'Latitude', 'Longitude', 'u', 'v', 'temperature', 'salinity'])<jupyter_text>Access salinity variable<jupyter_code>temp_c = f.get_variables_by_attributes(units='degC')[0]
sal = f.variables['salinity']
print(sal)
sal.dimensions, sal.shape<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>mt = f.variables['MT']
depth = f.variables['Depth']
x, y = f.variables['X'], f.variables['Y']
print(mt)
print(x)<jupyter_output><class 'netCDF4._netCDF4.Variable'>
float64 MT(MT)
    long_name: time
    units: days since 1900-12-31 00:00:00
    calendar: standard
    axis: T
unlimited dimensions: MT
current shape = (1,)
filling on, default _FillValue of 9.969209968386869e+36 used

<class 'netCDF4._netCDF4.Variable'>
int32 X(X)
    point_spacing: even
    axis: X
unlimited dimensions: 
current shape = (712,)
filling on, default _FillValue of -2147483647 used<jupyter_text>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.- 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).- Otherwise you get `IndexError: Index cannot be multidimensional` (?).  What region this dataset refers to? To plot it, we need the auxilary coordinate variables, `Latitude` and `Longitude` that define the 'real' data coordinates.<jupyter_code>lat = f.variables['Latitude']
lon = f.variables['Longitude']
print(lat)
# extract lat/lon values (in degrees) to numpy arrays
latvals = lat[:]; lonvals = lon[:] 
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
%matplotlib inline
fig, ax = plt.subplots(figsize=(12, 8),
                       subplot_kw=dict(projection=ccrs.PlateCarree()))
ax.coastlines('10m')
ax.set_title('Sea surface salinity', fontsize=24)
p = ax.contourf(lonvals, latvals, sal[0, 0, ...], cmap='viridis')
cb = plt.colorbar(p, orientation='horizontal', shrink=0.9, pad=0.05)

ax.set_extent([-180, -100, 40, 80])
ax.stock_img()

# Colorbar ticks fontsize
cb.ax.tick_params(labelsize=30)<jupyter_output><empty_output><jupyter_text>What is the sea surface temperature and salinity near St Lawrence island (63.202533N, 170.650635W)?- 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`<jupyter_code>def getclosest_ij(lats, lons, latpt, lonpt):
    """
    A function to find indices of the point closest
    (in squared distance) to return the given lat/lon value."""
    # 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)
pnt_lat, pnt_lon = (63.202533, -170.650635)
iy, ix = getclosest_ij(latvals, lonvals, pnt_lat, pnt_lon)
# Read values out of the netCDF file for temperature and salinity
print('{value:7.4f} {units}'.format(value=temp_c[0, 0, iy, ix], units=temp_c.units))
print('{value:7.4f} {units}'.format(value=sal[0, 0, iy, ix], units=sal.units))<jupyter_output>-1.1020 degC
32.0337 psu<jupyter_text>Simple multi-file aggregationWhat 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? Load the files<jupyter_code>files_wcard = os.path.join(os.path.pardir, 'data', 'uwnd*nc')
print(files_wcard)
mf = nc.MFDataset(files_wcard)<jupyter_output><empty_output><jupyter_text>Check time span<jupyter_code>times = mf.variables['time']
dates = nc.num2date(times[:], times.units)
print('starting date = {}'.format(dates[0]))
print('  ending date = {}'.format(dates[-1]))
u = mf.variables['uwnd']
print('u dimensions = {},\nu shape      = {}'.format(u.dimensions, u.shape))<jupyter_output>u dimensions = ('time', 'level', 'lat', 'lon'),
u shape      = (4018, 17, 73, 144)<jupyter_text>References This notebook is build upon the great materials of the [Unidata Python Workshop](https://unidata.github.io/unidata-python-workshop/):  * https://github.com/Unidata/unidata-python-workshop/blob/master/notebooks/netCDF-Reading.ipynb<jupyter_code>HTML(html)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1277749096541043, 'fraction_numerical': 0.034202374806401654, 'mean_word_length': 5.732406602953953, 'pattern_counts': {'":': 0, '<': 53, '<?xml version=': 0, '>': 53, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '324309', 'license': 'MIT', 'n_tokens_mistral': 2863, 'n_tokens_neox': 2686, 'n_words': 757, 'path': 'content/notebooks/2016-10-14-loading-netcdf-datasets.ipynb', 'repo_name': 'ueapy/ueapy.github.io', 'size': '7748'}