earthaccess 0.13.0
kerchunk 0.2.7
virtualizarr 1.3.2
xarray 2025.4.0
zarr 2.18.7
Note: you may need to restart the kernel to use updated packages.
# Built-in packagesimport osimport sys# Filesystem management import fsspecimport earthaccess# Data handlingimport xarray as xrimport numpy as np# Parallel computing import multiprocessingfrom dask import delayedimport dask.array as dafrom dask.distributed import Client# Otherimport matplotlib.pyplot as pltimport cartopy.crs as ccrsimport cartopy.feature as cfeature
Other Setup
xr.set_options( # display options for xarray objects display_expand_attrs=False, display_expand_coords=True, display_expand_data=True,)
<xarray.core.options.set_options at 0x7959127f1370>
1. Obtaining credentials and generating a VDS mapper
The first step is to find the S3 or https endpoints to the files. Handling access credentials to Earthdata and then finding the endpoints can be done a number of ways (e.g. using the requests, s3fs packages) but we use the earthaccess package for its ease of use.
# Get Earthdata credsearthaccess.login()
Enter your Earthdata Login username: edward.m.armstrong
Enter your Earthdata password: ········
<earthaccess.auth.Auth at 0x795a4da2d850>
In the next step we define a generic VDS wrapper and that can be wrapped into a function (no modification needed). The mapper will contain all the required access credentials and can be passed directly to xarray. In the future this step will likely be built into to earthaccess but for now we must define it in the notebook. The only inputs to the function are:
1. The link to the VDS reference file.
2. Whether or not you are accessing the data in the cloud in the same region as the data. For most beginning users this argument will be set to False.
def get_vds_mapper(vds_link, in_cloud_region=False):""" Produces a virtudal dataset mapper that can be passed to xarray. * vds_link: str, link to the mapper * in_cloud_region: bool, True if in cloud in the same region as the data, False otherwise. """if in_cloud_region: fs_data = earthaccess.get_s3_filesystem(daac="PODAAC") remote_protocol ="s3"else: fs_data = earthaccess.get_fsspec_https_session() remote_protocol ="https" storage_opts = {"fo": vds_link, "remote_protocol": remote_protocol, "remote_options": fs_data.storage_options} fs_ref = fsspec.filesystem('reference', **storage_opts)return fs_ref.get_mapper('')
2. Open reference files
OSTIA Daily Sea Surface Temperature (OSTIA-UKMO-L4-GLOB-REP-v2.0)
%%time# Open the OSTIA Reprocessed SST reference file (OSTIA-UKMO-L4-GLOB-REP-v2.0)vds_link ='https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-public/virtual_collections/OSTIA-UKMO-L4-GLOB-REP-v2.0/OSTIA-UKMO-L4-GLOB-REP-v2.0_virtual_https.json'vds_mapper = get_vds_mapper(vds_link, in_cloud_region=False)
CPU times: user 569 ms, sys: 166 ms, total: 735 ms
Wall time: 3.36 s
%%time## No modification needed. Read as xarray dataset.sst_ds = xr.open_dataset( vds_mapper, engine="zarr", chunks={}, backend_kwargs={"consolidated": False})print(sst_ds)sst_ds
<xarray.Dataset> Size: 11TB
Dimensions: (time: 15340, lat: 3600, lon: 7200)
Coordinates:
* lat (lat) float32 14kB -89.97 -89.93 -89.88 ... 89.93 89.97
* lon (lon) float32 29kB -180.0 -179.9 -179.9 ... 179.9 180.0
* time (time) datetime64[ns] 123kB 1982-01-01T12:00:00 ... 202...
Data variables:
analysed_sst (time, lat, lon) float64 3TB dask.array<chunksize=(1, 1200, 2400), meta=np.ndarray>
analysis_error (time, lat, lon) float64 3TB dask.array<chunksize=(1, 1200, 2400), meta=np.ndarray>
mask (time, lat, lon) float32 2TB dask.array<chunksize=(1, 1800, 3600), meta=np.ndarray>
sea_ice_fraction (time, lat, lon) float64 3TB dask.array<chunksize=(1, 1800, 3600), meta=np.ndarray>
Attributes: (39)
CPU times: user 600 ms, sys: 24.9 ms, total: 625 ms
Wall time: 644 ms
C.J. Donlon, M. Martin,J.D. Stark, J. Roberts-Jones, E. Fiedler, W. Wimmer. The operational sea surface temperature and sea ice analysis (OSTIA) system. Remote Sensing Environ., 116 (2012), pp. 140-158 http://dx.doi.org/10.1016/j.rse.2010.10.017
Donlon, C.J., Martin, M., Stark, J.D., Roberts-Jones, J., Fiedler, E., Wimmer, W., 2011. The Operational Sea Surface Temperature and Sea Ice Analysis (OSTIA). Remote Sensing of the Environment
institution :
UKMO
history :
Created from sst.nc; obs_anal.nc; seaice.nc
comment :
WARNING Some applications are unable to properly handle signed byte values. If values are encountered > 127, please subtract 256 from this reported value
license :
These data are available free of charge under the CMEMS data policy
RSS VAM 6-hour analyses using ERA-5 wind reanalyses as background. Satellite winds (vector and scalar) are assimilated using a variational analysis identical to that used in CCMP V1.1. The ERA5 winds are adjusted to match winds from scatteromters before use in the analysis.
institute_id :
RSS
institution :
Remote Sensing Systems (RSS)
source :
blend of satellite and ERA-5 winds
comment :
none
license :
This work is licensed under CC BY 4.0. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/
product_version :
3.1
standard_name_vocabulary :
NetCDF Climate and Forecast (CF) Metadata Convention
NASA Global Change Master Directory (GCMD) Science Keywords
platform_vocabulary :
GCMD platform keywords
instrument_vocabulary :
GCMD instrument keywords
references :
Mears, C.; Lee, T.; Ricciardulli, L.; Wang, X.; Wentz, F. Improving the Accuracy of the Cross-Calibrated Multi-Platform (CCMP) Ocean Vector Winds. Remote Sens. 2022, 14, 4230. https://doi.org/10.3390/rs14174230, Atlas, R., R. N. Hoffman, J. Ardizzone, S. M. Leidner, J. C. Jusem, D. K. Smith, D. Gombos, 2011: A cross-calibrated, multiplatform ocean surface wind velocity product for meteorological and oceanographic applications. Bull. Amer. Meteor. Soc., 92, 157-174. https://doi.org/10.1175/2010BAMS2946.1
# Determine and append the wind direction to the wind_ds dataset# Calculate wind direction in degrees using arctan2(-u, -v)wind_dir = np.arctan2(-wind_ds.uwnd, -wind_ds.vwnd) # radianswind_dir = np.degrees(wind_dir) # convert to degreeswind_dir = (wind_dir +360) %360# ensure values are between 0 and 360# Assign to datasetwind_ds['wind_direction'] = wind_dirwind_ds['wind_direction'].attrs.update({'units': 'degrees','standard_name': 'wind_from_direction','long_name': 'Wind Direction (meteorological)'})wind_ds
RSS VAM 6-hour analyses using ERA-5 wind reanalyses as background. Satellite winds (vector and scalar) are assimilated using a variational analysis identical to that used in CCMP V1.1. The ERA5 winds are adjusted to match winds from scatteromters before use in the analysis.
institute_id :
RSS
institution :
Remote Sensing Systems (RSS)
source :
blend of satellite and ERA-5 winds
comment :
none
license :
This work is licensed under CC BY 4.0. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/
product_version :
3.1
standard_name_vocabulary :
NetCDF Climate and Forecast (CF) Metadata Convention
NASA Global Change Master Directory (GCMD) Science Keywords
platform_vocabulary :
GCMD platform keywords
instrument_vocabulary :
GCMD instrument keywords
references :
Mears, C.; Lee, T.; Ricciardulli, L.; Wang, X.; Wentz, F. Improving the Accuracy of the Cross-Calibrated Multi-Platform (CCMP) Ocean Vector Winds. Remote Sens. 2022, 14, 4230. https://doi.org/10.3390/rs14174230, Atlas, R., R. N. Hoffman, J. Ardizzone, S. M. Leidner, J. C. Jusem, D. K. Smith, D. Gombos, 2011: A cross-calibrated, multiplatform ocean surface wind velocity product for meteorological and oceanographic applications. Bull. Amer. Meteor. Soc., 92, 157-174. https://doi.org/10.1175/2010BAMS2946.1
Region and time window selection
%%time# Define region # Gulf of Tehuantepec, MXlat_min =14.lat_max =15.lon_min =-96.lon_max =-95.lat_range = (lat_min, lat_max)lon_range = (lon_min, lon_max)# for datasets that use 0-360 deg lonlon_range_360 = (lon_min+360, lon_max+360)# Define the time slice start_date ='1993-01-01'#end_date = '1998-12-31'end_date ='2002-12-31'#end_date = '1994-12-31'time_range =(start_date, end_date)
CPU times: user 4 μs, sys: 0 ns, total: 4 μs
Wall time: 6.2 μs
Resample spatial SST and Wind Speed means to a common monthly time step and load into memory
%%time# A concise method to subset, resample, and average the xarray data all on one linesst_resample = sst_ds.analysed_sst.sel(lon=slice(*lon_range), lat=slice(*lat_range), time=slice(*time_range)).mean(["lat", "lon"]).resample(time="1ME").mean().load()wind_speed_resample = wind_ds.ws.sel(longitude=slice(*lon_range_360), latitude=slice(*lat_range), time=slice(*time_range)).mean(["latitude", "longitude"]).resample(time="1ME").mean().load()
CPU times: user 5min 30s, sys: 33.9 s, total: 6min 4s
Wall time: 3min 38s
Resample the wind directions means to the same common monthly time step and find the upwelling favorable wind direction and speed times
%%time# Resample u and v winds (monthly mean)u_resample = wind_ds.uwnd.sel( longitude=slice(*lon_range_360), latitude=slice(*lat_range), time=slice(*time_range)).mean(["latitude", "longitude"]).resample(time="1ME").mean().load()v_resample = wind_ds.vwnd.sel( longitude=slice(*lon_range_360), latitude=slice(*lat_range), time=slice(*time_range)).mean(["latitude", "longitude"]).resample(time="1ME").mean().load()# Resample wind directionwind_dir_resample = wind_ds.wind_direction.sel( longitude=slice(*lon_range_360), latitude=slice(*lat_range), time=slice(*time_range)).mean(["latitude", "longitude"]).resample(time="1ME").mean().load()
CPU times: user 23min 49s, sys: 1min 2s, total: 24min 51s
Wall time: 8min 52s
# Identify upwelling favorable wind direction times (0–120 degrees) and speeds of at least 6.8 m/sfavorable_mask = (wind_dir_resample >=0 ) & (wind_dir_resample <=120) & (wind_speed_resample >=6.8)favorable_times = wind_dir_resample.time[favorable_mask]
Check with a plot
fig, ax1 = plt.subplots(figsize=(14, 6))# SST Plotcolor ='tab:blue'ax1.set_xlabel('Year')ax1.set_ylabel('SST (°C)', color=color, fontsize='x-large')ax1.set_xlabel('Year', fontsize='x-large')ax1.plot(sst_resample['time'], (sst_resample -273.15), linewidth=2, color=color, label='SST')ax1.tick_params(axis='y', labelcolor=color, size=8, labelsize=12)ax1.tick_params(axis='x', size=8, labelsize=12)# Wind Speed Plotax2 = ax1.twinx()color ='tab:red'ax2.set_ylabel('Wind Speed (m/s)', color=color, fontsize='x-large')ax2.plot(wind_speed_resample['time'], wind_speed_resample, color=color, label='Wind Speed')ax2.tick_params(axis='y', labelcolor=color, size=8, labelsize=12)# Add vertical shaded regions for favorable wind direction periodsfor ft in favorable_times.values: ax1.axvspan(ft - np.timedelta64(15, 'D'), ft + np.timedelta64(15, 'D'), color='lightgray', alpha=0.5)# Add wind barbs in lower panel (directional vectors)# Normalize u and v components for visual lengthbarb_skip =2# Plot every nth point to reduce clutter if neededbarb_scale =30# Adjust scale for clarity# Only plot barbs every few monthsbarb_time = u_resample.time.values[::barb_skip]u_plot = u_resample.values[::barb_skip]v_plot = v_resample.values[::barb_skip]# Set position for barbs slightly below wind speed linebarb_y = np.full_like(u_plot, wind_speed_resample.min().item() -1)# Barbs on ax2 (wind axis)ax2.barbs(barb_time, barb_y, u_plot, v_plot, length=6, pivot='middle', barb_increments=dict(half=1, full=2, flag=5), color='k', linewidth=0.5)# Add legendlines_1, labels_1 = ax1.get_legend_handles_labels()lines_2, labels_2 = ax2.get_legend_handles_labels()ax2.legend(lines_1 + lines_2, labels_1 + labels_2, loc='upper left')plt.title("SST, Ocean Surface Wind Speed and Direction with Upwelling Favorable Periods Shaded", fontsize='x-large')fig.tight_layout()plt.show()
Zonal (east-west) component of ocean velocity. Note: EVEL is calculated by interpolating the model's x and y components of ocean velocity (UVEL and VVEL)to tracer cell centers and then finding the zonal component of the interpolated vectors. It is not recommended to use EVEL and NVEL for volume budget calculations because interpolating UVEL and VVEL from the model grid to the lat-lon grid introduces errors. Perform volume budget calculations with UVELMASS and VVELMASS on the native model grid.
Meridional (north-south) component of ocean velocity. Note: NVEL is calculated by interpolating the model's x and y components of ocean velocity (UVEL and VVEL) to tracer cell centers and then finding the meridional component of the interpolated vectors. It is not recommended to use EVEL and NVEL for volume budget calculations because interpolating UVEL and VVEL from the model grid to the lat-lon grid introduces errors. Perform volume budget calculations with UVELMASS and VVELMASS on the native model grid.
Vertical velocity in the +z direction at the top face of the grid cell. Note: in the Arakawa-C grid used in ECCO V4r4, vertical velocities are staggered relative to the tracer cell centers with values at the TOP and BOTTOM faces of each grid cell.
This research was carried out by the Jet Propulsion Laboratory, managed by the California Institute of Technology under a contract with the National Aeronautics and Space Administration.
author :
Ian Fenty and Ou Wang
cdm_data_type :
Grid
comment :
Fields provided on a regular lat-lon grid. They have been mapped to the regular lat-lon grid from the original ECCO lat-lon-cap 90 (llc90) native model grid.
Conventions :
CF-1.8, ACDD-1.3
coordinates_comment :
Note: the global 'coordinates' attribute describes auxillary coordinates.
creator_email :
ecco-group@mit.edu
creator_institution :
NASA Jet Propulsion Laboratory (JPL)
creator_name :
ECCO Consortium
creator_type :
group
creator_url :
https://ecco-group.org
geospatial_bounds_crs :
EPSG:4326
geospatial_lat_max :
90.0
geospatial_lat_min :
-90.0
geospatial_lat_resolution :
0.5
geospatial_lat_units :
degrees_north
geospatial_lon_max :
180.0
geospatial_lon_min :
-180.0
geospatial_lon_resolution :
0.5
geospatial_lon_units :
degrees_east
geospatial_vertical_max :
0.0
geospatial_vertical_min :
-6134.5
geospatial_vertical_positive :
up
geospatial_vertical_resolution :
variable
geospatial_vertical_units :
meter
history :
Inaugural release of an ECCO Central Estimate solution to PO.DAAC
id :
10.5067/ECG5D-OVE44
institution :
NASA Jet Propulsion Laboratory (JPL)
instrument_vocabulary :
GCMD instrument keywords
keywords_vocabulary :
NASA Global Change Master Directory (GCMD) Science Keywords
NASA Physical Oceanography, Cryosphere, Modeling, Analysis, and Prediction (MAP)
project :
Estimating the Circulation and Climate of the Ocean (ECCO)
publisher_email :
podaac@podaac.jpl.nasa.gov
publisher_institution :
PO.DAAC
publisher_name :
Physical Oceanography Distributed Active Archive Center (PO.DAAC)
publisher_type :
institution
publisher_url :
https://podaac.jpl.nasa.gov
references :
ECCO Consortium, Fukumori, I., Wang, O., Fenty, I., Forget, G., Heimbach, P., & Ponte, R. M. 2020. Synopsis of the ECCO Central Production Global Ocean and Sea-Ice State Estimate (Version 4 Release 4). doi:10.5281/zenodo.3765928
source :
The ECCO V4r4 state estimate was produced by fitting a free-running solution of the MITgcm (checkpoint 66g) to satellite and in situ observational data in a least squares sense using the adjoint method
standard_name_vocabulary :
NetCDF Climate and Forecast (CF) Metadata Convention
summary :
This dataset provides daily-averaged ocean velocity interpolated to a regular 0.5-degree grid from the ECCO Version 4 Release 4 (V4r4) ocean and sea-ice state estimate. Estimating the Circulation and Climate of the Ocean (ECCO) state estimates are dynamically and kinematically-consistent reconstructions of the three-dimensional, time-evolving ocean, sea-ice, and surface atmospheric states. ECCO V4r4 is a free-running solution of a global, nominally 1-degree configuration of the MIT general circulation model (MITgcm) that has been fit to observations in a least-squares sense. Observational data constraints used in V4r4 include sea surface height (SSH) from satellite altimeters [ERS-1/2, TOPEX/Poseidon, GFO, ENVISAT, Jason-1,2,3, CryoSat-2, and SARAL/AltiKa]; sea surface temperature (SST) from satellite radiometers [AVHRR], sea surface salinity (SSS) from the Aquarius satellite radiometer/scatterometer, ocean bottom pressure (OBP) from the GRACE satellite gravimeter; sea-ice concentration from satellite radiometers [SSM/I and SSMIS], and in-situ ocean temperature and salinity measured with conductivity-temperature-depth (CTD) sensors and expendable bathythermographs (XBTs) from several programs [e.g., WOCE, GO-SHIP, Argo, and others] and platforms [e.g., research vessels, gliders, moorings, ice-tethered profilers, and instrumented pinnipeds]. V4r4 covers the period 1992-01-01T12:00:00 to 2018-01-01T00:00:00.
Resample the ocean vertical velocity means to the same common monthly time step
%%time# Do the vertical velocity resampling for the first depth level (5 m below the surface)depth_min =-5depth_max =0depth_range = (depth_min, depth_max)vertical_vel_resample = ecco_current_ds.WVEL.sel(longitude=slice(*lon_range), latitude=slice(*lat_range), time=slice(*time_range), Z=slice(0,-5)).mean(["latitude", "longitude"]).resample(time="1ME").mean().load()
CPU times: user 1min 35s, sys: 3.59 s, total: 1min 38s
Wall time: 47.4 s
fig, ax1 = plt.subplots(figsize=(14, 6))# SST Plotcolor ='tab:blue'ax1.set_xlabel('Year')ax1.set_ylabel('SST (°C)', color=color, fontsize='x-large')ax1.set_xlabel('Year', fontsize='x-large')ax1.plot(sst_resample['time'], (sst_resample -273.15), linewidth=2, color=color, label='SST')ax1.tick_params(axis='y', labelcolor=color, size=8, labelsize=12)ax1.tick_params(axis='x', size=8, labelsize=12)# Ocean Vertical Velocity Plotax2 = ax1.twinx()color ='tab:purple'ax2.set_ylabel('Vertical Velocity (m/s)', color=color, fontsize='x-large')ax2.plot(vertical_vel_resample['time'], vertical_vel_resample, color=color, label='Vertical Velocity (5m)')ax2.tick_params(axis='y', labelcolor=color, size=8, labelsize=12)# Mask for the positive (upwelling) vertical velocitiespositive_mask = (vertical_vel_resample >0)# Plot positive velocities (in a highlighted style)ax2.plot(vertical_vel_resample['time'].where(positive_mask), vertical_vel_resample.where(positive_mask), color=color, linewidth=3, label='Positive Vertical Velocity (5m)')# Add vertical shaded regions for favorable wind direction periods from previous plotfor ft in favorable_times.values: ax1.axvspan(ft - np.timedelta64(15, 'D'), ft + np.timedelta64(15, 'D'), color='lightgray', alpha=0.5)# Add legendlines_1, labels_1 = ax1.get_legend_handles_labels()lines_2, labels_2 = ax2.get_legend_handles_labels()ax2.legend(lines_1 + lines_2, labels_1 + labels_2, loc='upper left')plt.title("SST and Ocean Vertical Velocity with Upwelling Favorable Periods Shaded", fontsize='x-large')fig.tight_layout()plt.show()
ECCO Ocean Temperature and Salinity - Daily Mean 0.5 Degree (Version 4 Release 4) dataset (ECCO_L4_TEMP_SALINITY_05DEG_DAILY_V4R4)
%%time# Open the ECCO Ocean Temperature and Salinity - Daily Mean 0.5 Degree (Version 4 Release 4) dataset (ECCO_L4_TEMP_SALINITY_05DEG_DAILY_V4R4)vds_link ='https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-public/virtual_collections/ECCO_L4_TEMP_SALINITY_05DEG_DAILY_V4R4/ECCO_L4_TEMP_SALINITY_05DEG_DAILY_V4R4_virtual_https.json'vds_mapper = get_vds_mapper(vds_link, in_cloud_region=False)ecco_temp_salt_ds = xr.open_dataset( vds_mapper, engine="zarr", chunks={}, backend_kwargs={"consolidated": False})print(ecco_temp_salt_ds)ecco_temp_salt_ds
Defined using CF convention 'Sea water salinity is the salt content of sea water, often on the Practical Salinity Scale of 1978. However, the unqualified term 'salinity' is generic and does not necessarily imply any particular method of calculation. The units of salinity are dimensionless and the units attribute should normally be given as 1e-3 or 0.001 i.e. parts per thousand.' see https://cfconventions.org/Data/cf-standard-names/73/build/cf-standard-name-table.html
Sea water potential temperature is the temperature a parcel of sea water would have if moved adiabatically to sea level pressure. Note: the equation of state is a modified UNESCO formula by Jackett and McDougall (1995), which uses the model variable potential temperature as input assuming a horizontally and temporally constant pressure of $p_0=-g ho_{0} z$.
This research was carried out by the Jet Propulsion Laboratory, managed by the California Institute of Technology under a contract with the National Aeronautics and Space Administration.
author :
Ian Fenty and Ou Wang
cdm_data_type :
Grid
comment :
Fields provided on a regular lat-lon grid. They have been mapped to the regular lat-lon grid from the original ECCO lat-lon-cap 90 (llc90) native model grid.
Conventions :
CF-1.8, ACDD-1.3
coordinates_comment :
Note: the global 'coordinates' attribute describes auxillary coordinates.
creator_email :
ecco-group@mit.edu
creator_institution :
NASA Jet Propulsion Laboratory (JPL)
creator_name :
ECCO Consortium
creator_type :
group
creator_url :
https://ecco-group.org
geospatial_bounds_crs :
EPSG:4326
geospatial_lat_max :
90.0
geospatial_lat_min :
-90.0
geospatial_lat_resolution :
0.5
geospatial_lat_units :
degrees_north
geospatial_lon_max :
180.0
geospatial_lon_min :
-180.0
geospatial_lon_resolution :
0.5
geospatial_lon_units :
degrees_east
geospatial_vertical_max :
0.0
geospatial_vertical_min :
-6134.5
geospatial_vertical_positive :
up
geospatial_vertical_resolution :
variable
geospatial_vertical_units :
meter
history :
Inaugural release of an ECCO Central Estimate solution to PO.DAAC
id :
10.5067/ECG5D-OTS44
institution :
NASA Jet Propulsion Laboratory (JPL)
instrument_vocabulary :
GCMD instrument keywords
keywords_vocabulary :
NASA Global Change Master Directory (GCMD) Science Keywords
NASA Physical Oceanography, Cryosphere, Modeling, Analysis, and Prediction (MAP)
project :
Estimating the Circulation and Climate of the Ocean (ECCO)
publisher_email :
podaac@podaac.jpl.nasa.gov
publisher_institution :
PO.DAAC
publisher_name :
Physical Oceanography Distributed Active Archive Center (PO.DAAC)
publisher_type :
institution
publisher_url :
https://podaac.jpl.nasa.gov
references :
ECCO Consortium, Fukumori, I., Wang, O., Fenty, I., Forget, G., Heimbach, P., & Ponte, R. M. 2020. Synopsis of the ECCO Central Production Global Ocean and Sea-Ice State Estimate (Version 4 Release 4). doi:10.5281/zenodo.3765928
source :
The ECCO V4r4 state estimate was produced by fitting a free-running solution of the MITgcm (checkpoint 66g) to satellite and in situ observational data in a least squares sense using the adjoint method
standard_name_vocabulary :
NetCDF Climate and Forecast (CF) Metadata Convention
summary :
This dataset provides daily-averaged ocean potential temperature and salinity interpolated to a regular 0.5-degree grid from the ECCO Version 4 Release 4 (V4r4) ocean and sea-ice state estimate. Estimating the Circulation and Climate of the Ocean (ECCO) state estimates are dynamically and kinematically-consistent reconstructions of the three-dimensional, time-evolving ocean, sea-ice, and surface atmospheric states. ECCO V4r4 is a free-running solution of a global, nominally 1-degree configuration of the MIT general circulation model (MITgcm) that has been fit to observations in a least-squares sense. Observational data constraints used in V4r4 include sea surface height (SSH) from satellite altimeters [ERS-1/2, TOPEX/Poseidon, GFO, ENVISAT, Jason-1,2,3, CryoSat-2, and SARAL/AltiKa]; sea surface temperature (SST) from satellite radiometers [AVHRR], sea surface salinity (SSS) from the Aquarius satellite radiometer/scatterometer, ocean bottom pressure (OBP) from the GRACE satellite gravimeter; sea-ice concentration from satellite radiometers [SSM/I and SSMIS], and in-situ ocean temperature and salinity measured with conductivity-temperature-depth (CTD) sensors and expendable bathythermographs (XBTs) from several programs [e.g., WOCE, GO-SHIP, Argo, and others] and platforms [e.g., research vessels, gliders, moorings, ice-tethered profilers, and instrumented pinnipeds]. V4r4 covers the period 1992-01-01T12:00:00 to 2018-01-01T00:00:00.
time_coverage_duration :
P1D
time_coverage_resolution :
P1D
title :
ECCO Ocean Temperature and Salinity - Daily Mean 0.5 Degree (Version 4 Release 4)
Resample the ocean surface salinity means to the same common monthly time step
%%time# Do the ocean salinity resampling for the first depth level (5 m below the surface)depth_min =-5depth_max =0depth_range = (depth_min, depth_max)salt_resample = ecco_temp_salt_ds.SALT.sel(longitude=slice(*lon_range), latitude=slice(*lat_range), time=slice(*time_range), Z=slice(0,-5)).mean(["latitude", "longitude"]).resample(time="1ME").mean().load()
CPU times: user 1min 20s, sys: 2.78 s, total: 1min 22s
Wall time: 44.9 s
fig, ax1 = plt.subplots(figsize=(14, 6))# SST Plotcolor ='tab:blue'ax1.set_xlabel('Year')ax1.set_ylabel('SST (°C)', color=color, fontsize='x-large')ax1.set_xlabel('Year', fontsize='x-large')ax1.plot(sst_resample['time'], (sst_resample -273.15), linewidth=2, color=color, label='SST')ax1.tick_params(axis='y', labelcolor=color, size=8, labelsize=12)ax1.tick_params(axis='x', size=8, labelsize=12)# Ocean Surface Salinity Plotax2 = ax1.twinx()color ='tab:orange'ax2.set_ylabel('SSS (ppt)', color=color, fontsize='x-large')ax2.plot(salt_resample['time'], salt_resample, color=color, label='Sea Surface Salinity (5m)')ax2.tick_params(axis='y', labelcolor=color, size=8, labelsize=12)# Add vertical shaded regions for favorable wind direction periods from previous plotfor ft in favorable_times.values: ax1.axvspan(ft - np.timedelta64(15, 'D'), ft + np.timedelta64(15, 'D'), color='lightgray', alpha=0.5)# Add legendlines_1, labels_1 = ax1.get_legend_handles_labels()lines_2, labels_2 = ax2.get_legend_handles_labels()ax2.legend(lines_1 + lines_2, labels_1 + labels_2, loc='upper left')plt.title("SST and Sea Surface Salinity (SSS) with Upwelling Favorable Periods Shaded", fontsize='x-large')fig.tight_layout()plt.show()
Use parallel processing to speed the resample() implementations
# Set up parallel processing# Check how many cpu's are on this VM:print("CPU count =", multiprocessing.cpu_count())# Start up cluster and print some information about it:client = Client(n_workers=32, threads_per_worker=1)print(client.cluster)print("View any work being done on the cluster here", client.dashboard_link)
CPU count = 32
LocalCluster(c6deaf14, 'tcp://127.0.0.1:35297', workers=32, threads=32, memory=122.27 GiB)
View any work being done on the cluster here https://cluster-nwbxl.dask.host/jupyter/proxy/8787/status
%%time# Execute Tasks: Convert to weekly (1W) monthly (1ME) or annual means (1YE). # use load() to load into memory the result# A concise method to subset, resample, and average the xarray data all on one linesst_resample = sst_ds.analysed_sst.sel(lon=slice(*lon_range), lat=slice(*lat_range), time=slice(*time_range)).mean(["lat", "lon"]).resample(time="1ME").mean().load()wind_speed_resample = wind_ds.ws.sel(longitude=slice(*lon_range_360), latitude=slice(*lat_range), time=slice(*time_range)).mean(["latitude", "longitude"]).resample(time="1ME").mean().load()u_resample = wind_ds.uwnd.sel( longitude=slice(*lon_range_360), latitude=slice(*lat_range), time=slice(*time_range)).mean(["latitude", "longitude"]).resample(time="1ME").mean().load()v_resample = wind_ds.vwnd.sel( longitude=slice(*lon_range_360), latitude=slice(*lat_range), time=slice(*time_range)).mean(["latitude", "longitude"]).resample(time="1ME").mean().load()wind_dir_resample = wind_ds.wind_direction.sel( longitude=slice(*lon_range_360), latitude=slice(*lat_range), time=slice(*time_range)).mean(["latitude", "longitude"]).resample(time="1ME").mean().load()# vertical velocity# surface salinity# Trigger computation with parallel execution....optionalsst_data, wind_speed_data, u_data, v_data, wind_dir_data = da.compute(sst_resample, wind_speed_resample, u_resample, v_resample, wind_dir_resample)
CPU times: user 4min 39s, sys: 47.4 s, total: 5min 26s
Wall time: 13min 50s
Check with a plot
# Identify favorable wind direction times (0–120 degrees) and speeds of at least 6.8 m/sfavorable_mask = (wind_dir_resample >=0 ) & (wind_dir_resample <=120) & (wind_speed_resample >=6.8)favorable_times = wind_dir_resample.time[favorable_mask]
fig, ax1 = plt.subplots(figsize=(14, 6))# SST Plotcolor ='tab:blue'ax1.set_xlabel('Year')ax1.set_ylabel('SST (°K)', color=color)ax1.plot(sst_resample['time'], sst_resample, linewidth=2, color=color, label='SST')ax1.tick_params(axis='y', labelcolor=color)# Wind Speed Plotax2 = ax1.twinx()color ='tab:red'ax2.set_ylabel('Wind Speed (m/s)', color=color)ax2.plot(wind_speed_resample['time'], wind_speed_resample, color=color, label='Wind Speed')ax2.tick_params(axis='y', labelcolor=color)# Add vertical shaded regions for favorable wind direction periodsfor ft in favorable_times.values: ax1.axvspan(ft - np.timedelta64(15, 'D'), ft + np.timedelta64(15, 'D'), color='lightgray', alpha=0.3)# Add wind barbs in lower panel (directional vectors)# Normalize u and v components for visual lengthbarb_skip =2# Plot every nth point to reduce clutter if neededbarb_scale =30# Adjust scale for clarity# Only plot barbs every few monthsbarb_time = u_resample.time.values[::barb_skip]u_plot = u_resample.values[::barb_skip]v_plot = v_resample.values[::barb_skip]# Set position for barbs slightly below wind speed linebarb_y = np.full_like(u_plot, wind_speed_resample.min().item() -1)# Barbs on ax2 (wind axis)ax2.barbs(barb_time, barb_y, u_plot, v_plot, length=5, pivot='middle', barb_increments=dict(half=1, full=2, flag=5), color='k', linewidth=0.5)# Add legendlines_1, labels_1 = ax1.get_legend_handles_labels()lines_2, labels_2 = ax2.get_legend_handles_labels()ax2.legend(lines_1 + lines_2, labels_1 + labels_2, loc='upper left')plt.title("SST, Ocean Surface Wind Speed, and Wind Direction (with Favorable Periods Shaded)")fig.tight_layout()plt.show()
