<jupyter_start><jupyter_text>Direct access to ECCO data in S3 (from *us-west-2*) Getting StartedIn this notebook, we will access monthly sea surface height from ECCO V4r4 ([10.5067/ECG5D-SSH44](https://doi.org/10.5067/ECG5D-SSH44)). The data are provided as a time series of monthly netCDFs on a 0.5-degree latitude/longitude grid. We will access the data from inside the AWS cloud (*us-west-2* region, specifically) and load a time series made of multiple netCDF datasets into a single xarray dataset. This approach leverages S3 native protocols for efficient access to the data. Requirements AWSThis notebook should be running in an EC2 instance in AWS region *us-west-2*, as previously mentioned. We recommend using an EC2 with at least 8GB of memory available. The notebook was developed and tested using a *t2.large* instance (2 cpus; 8GB memory). Python 3Most of these imports are from the Python standard library. However, you will need to install these packages into your Python 3 environment if you have not already done so:* s3fs* requests* pandas* xarray* matplotlib* cartopy<jupyter_code>import s3fs
import requests
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeat
from json import dumps
from io import StringIO
from os.path import dirname, join
from IPython.display import HTML

plt.rcParams.update({'font.size': 14})<jupyter_output><empty_output><jupyter_text>Make a folder to write some outputs, if needed:<jupyter_code>!mkdir -p outputs/<jupyter_output><empty_output><jupyter_text>InputsConfigure one input: the *ShortName* of the desired dataset from ECCO V4r4. In this case it's the following string that unique identifies the collection of monthly, 0.5-degree sea surface height data.<jupyter_code>ShortName = "ECCO_L4_SSH_05DEG_MONTHLY_V4R4"<jupyter_output><empty_output><jupyter_text>Earthdata LoginYou should have a *.netrc* file set up like:```machine urs.earthdata.nasa.gov login  password ``` Direct access from S3Set up an s3fs session for authneticated access to ECCO netCDF files in s3:<jupyter_code>def begin_s3_direct_access(url: str="https://archive.podaac.earthdata.nasa.gov/s3credentials"):
    response = requests.get(url).json()
    return s3fs.S3FileSystem(key=response['accessKeyId'],
                             secret=response['secretAccessKey'],
                             token=response['sessionToken'],
                             client_kwargs={'region_name':'us-west-2'})

fs = begin_s3_direct_access()

type(fs)<jupyter_output><empty_output><jupyter_text>Datasets sea surface height (0.5-degree gridded, monthly)[ECCO_L4_SSH_05DEG_MONTHLY_V4R4](https://podaac.jpl.nasa.gov/dataset/ECCO_L4_SSH_05DEG_MONTHLY_V4R4)Get a list of netCDF files located at the S3 path corresponding to the ECCO V4r4 monthly sea surface height dataset on the 0.5-degree latitude/longitude grid.<jupyter_code>ssh_Files = fs.glob(join("podaac-ops-cumulus-protected/", ShortName, "*2015*.nc"))

len(ssh_Files)<jupyter_output><empty_output><jupyter_text>Open with the netCDF files using the *s3fs* package, then load them all at once into a concatenated *xarray* dataset.<jupyter_code>ssh_Dataset = xr.open_mfdataset(
    paths=[fs.open(f) for f in ssh_Files],
    combine='by_coords',
    mask_and_scale=True,
    decode_cf=True,
    chunks={'latitude': 60,   # These were chosen arbitrarily. You must specify 
            'longitude': 120, # chunking that is suitable to the data and target
            'time': 100}      # analysis.
)

ssh = ssh_Dataset.SSH

print(ssh)<jupyter_output><xarray.DataArray 'SSH' (time: 12, latitude: 360, longitude: 720)>
dask.array<concatenate, shape=(12, 360, 720), dtype=float32, chunksize=(1, 60, 120), chunktype=numpy.ndarray>
Coordinates:
  * time       (time) datetime64[ns] 2015-01-16T12:00:00 ... 2015-12-16T12:00:00
  * latitude   (latitude) float32 -89.75 -89.25 -88.75 ... 88.75 89.25 89.75
  * longitude  (longitude) float32 -179.8 -179.2 -178.8 ... 178.8 179.2 179.8
Attributes:
    coverage_content_type:  modelResult
    long_name:              Dynamic sea surface height anomaly
    standard_name:          sea_surface_height_above_geoid
    units:                  m
    comment:                Dynamic sea surface height anomaly above the geoi...
    valid_min:              [-1.88057721]
    valid_max:              [1.42077196]<jupyter_text>Plot the gridded *sea surface height* time seriesBut only the timesteps beginning in 2015:<jupyter_code>ssh_after_201x = ssh[ssh['time.year']>=2015,:,:]

print(ssh_after_201x)<jupyter_output><xarray.DataArray 'SSH' (time: 12, latitude: 360, longitude: 720)>
dask.array<concatenate, shape=(12, 360, 720), dtype=float32, chunksize=(1, 60, 120), chunktype=numpy.ndarray>
Coordinates:
  * time       (time) datetime64[ns] 2015-01-16T12:00:00 ... 2015-12-16T12:00:00
  * latitude   (latitude) float32 -89.75 -89.25 -88.75 ... 88.75 89.25 89.75
  * longitude  (longitude) float32 -179.8 -179.2 -178.8 ... 178.8 179.2 179.8
Attributes:
    coverage_content_type:  modelResult
    long_name:              Dynamic sea surface height anomaly
    standard_name:          sea_surface_height_above_geoid
    units:                  m
    comment:                Dynamic sea surface height anomaly above the geoi...
    valid_min:              [-1.88057721]
    valid_max:              [1.42077196]<jupyter_text>Plot the grid for the first time step using a Robinson projection. Define a helper function for consistency throughout the notebook:<jupyter_code>def make_figure(proj):
    fig = plt.figure(figsize=(16,6))
    ax = fig.add_subplot(1, 1, 1, projection=proj)
    ax.add_feature(cfeat.LAND)
    ax.add_feature(cfeat.OCEAN)
    ax.add_feature(cfeat.COASTLINE)
    ax.add_feature(cfeat.BORDERS, linestyle='dotted')
    return fig, ax

fig, ax = make_figure(proj=ccrs.Robinson())

ssh_after_201x.isel(time=0).plot(ax=ax, transform=ccrs.PlateCarree(), cmap='Spectral_r')<jupyter_output><empty_output><jupyter_text>Now plot the whole time series (post-2010) in an animation and write it to an mp4 file called *ecco_monthly_ssh_grid_2015_to_x.mp4*:<jupyter_code>def get_animation(var, cmap: str="Spectral_r"):
    """Get time series animation for input xarray dataset"""

    def draw_map(i: int, add_colorbar: bool):
        data = var[i]
        m = data.plot(ax=ax, 
                      transform=ccrs.PlateCarree(),
                      add_colorbar=add_colorbar,
                      vmin=var.valid_min, 
                      vmax=var.valid_max,
                      cmap=cmap)
        plt.title(str(data.time.values)[:7])
        return m

    def init():
        return draw_map(0, add_colorbar=True)
    
    def animate(i):
        return draw_map(i, add_colorbar=False)

    return init, animate<jupyter_output><empty_output><jupyter_text>Now make the animation using the function:<jupyter_code>fig, ax = make_figure(proj=ccrs.Robinson())

init, animate = get_animation(ssh_after_201x)

ani = animation.FuncAnimation(fig=fig, 
                              func=animate, 
                              frames=ssh_after_201x.time.size, 
                              init_func=init, 
                              interval=0.2, 
                              blit=False, 
                              repeat=False)

# Now save the animation to an MP4 file:
ani.save('outputs/ecco_monthly_ssh_grid_2015_to_x.mp4', writer=animation.FFMpegWriter(fps=8))

plt.close(fig)<jupyter_output><empty_output><jupyter_text>Render the animation in the ipynb:<jupyter_code>#HTML(ani.to_html5_video())<jupyter_output><empty_output><jupyter_text>tflux (0.5-degree gridded, daily)Now we will do something similar to access daily, gridded (0.5-degree) ocean and sea-ice surface heat fluxes ([10.5067/ECG5D-HEA44](https://doi.org/10.5067/ECG5D-HEA44)). Read more about the dataset and the rest of the ECCO V4r4 product suite on the [PO.DAAC Web Portal](https://podaac.jpl.nassa.gov). Use a "glob" pattern when listing the S3 bucket contents such that only netCDFs from January 2015 are represented in the resulting list of paths.<jupyter_code>tflux_Files = fs.glob(join("podaac-ops-cumulus-protected/", "ECCO_L4_HEAT_FLUX_05DEG_DAILY_V4R4", "*2015-01*.nc"))

len(tflux_Files)<jupyter_output><empty_output><jupyter_text>Now open them all as one xarray dataset just like before. Open and pass the 365 netCDF files to the *xarray.open_mfdataset* constructor so that we can operate on them as a single aggregated dataset.<jupyter_code>tflux_Dataset = xr.open_mfdataset(
    paths=[fs.open(f) for f in tflux_Files],
    combine='by_coords',
    mask_and_scale=True,
    decode_cf=True,
    chunks={'latitude': 60,   # These were chosen arbitrarily. You must specify 
            'longitude': 120, # chunking that is suitable to the data and target
            'time': 100}      # analysis.
)

tflux = tflux_Dataset.TFLUX

print(tflux)<jupyter_output><xarray.DataArray 'TFLUX' (time: 31, latitude: 360, longitude: 720)>
dask.array<concatenate, shape=(31, 360, 720), dtype=float32, chunksize=(1, 60, 120), chunktype=numpy.ndarray>
Coordinates:
  * time       (time) datetime64[ns] 2015-01-01T12:00:00 ... 2015-01-31T12:00:00
  * latitude   (latitude) float32 -89.75 -89.25 -88.75 ... 88.75 89.25 89.75
  * longitude  (longitude) float32 -179.8 -179.2 -178.8 ... 178.8 179.2 179.8
Attributes:
    coverage_content_type:  modelResult
    direction:              >0 increases potential temperature (THETA)
    long_name:              Rate of change of ocean heat content per m2 accou...
    units:                  W m-2
    comment:                The rate of change of ocean heat content due to h...
    valid_min:              [-1713.51220703]
    valid_max:              [870.31304932]<jupyter_text>Select a region over the Gulf of Mexico and spatially subset it from the larger dataset by slicing on the *latitude* and *longitude* axes.<jupyter_code>tflux_gom = tflux.sel(latitude=slice(15, 40), 
                      longitude=slice(-105, -70))

print(tflux_gom.shape)
tflux_gom.isel(time=0).plot()<jupyter_output><empty_output><jupyter_text>Plot the Jan 2015 surface heat flux as a gridded time series animation over the GOM study region.<jupyter_code>fig, ax = make_figure(proj=ccrs.Mercator())

ax.coastlines()
ax.set_extent([tflux_gom.longitude.min(), 
               tflux_gom.longitude.max(), 
               tflux_gom.latitude.min(), 
               tflux_gom.latitude.max()])

init, animate = get_animation(tflux_gom, cmap="RdBu")

# Plot a time series animation write it to an mp4 file:
ani = animation.FuncAnimation(fig=fig, 
                              func=animate, 
                              frames=tflux_gom.time.size, 
                              init_func=init, 
                              interval=0.2, 
                              blit=False, 
                              repeat=False)

ani.save('outputs/ecco_daily_tflux_gom_2015.mp4', writer=animation.FFMpegWriter(fps=8))

plt.close(fig)<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup