1. How to work with satellite data in Matlab

This tutorial will show the steps to grab data in ERDDAP from Matlab, how to work with NetCDF files in Matlab, and how to make some maps and time-series of sea surface temperature (SST) around the main Hawaiian islands.
Note: The mapping code in this tutorial uses the Mapping Toolbox. If you don't have access to the Mapping Toolbox, you can still create things like filled contour plots of your data using 'contourf' or an open source package like M_map (https://www.eoas.ubc.ca/~rich/map.html) or the Climate Data Toolbox for Matlab (https://github.com/chadagreene/CDT).

Downloading data in Matlab

Because ERDDAP includes RESTful services, you can download data listed on any ERDDAP platform from Matlab using the URL structure. For example, the following page allows you to subset monthly SST data: https://oceanwatch.pifsc.noaa.gov/erddap/griddap/CRW_sst_v3_1_monthly.html
Select your region and date range of interest, then select the '.nc' (NetCDF) file type and click on "Just Generate the URL".
In this specific example, the URL we generated is:
https://oceanwatch.pifsc.noaa.gov/erddap/griddap/CRW_sst_v3_1_monthly.nc?sea_surface_temperature%5B(2018-1-31T12:00:00Z):1:(2018-12-31T12:00:00Z)%5D%5B(17):1:(30)%5D%5B(195):1:(210)%5D
You can also edit this URL manually.

Importing the data in Matlab

In Matlab, run the following code to view details about the data using a portion of the generated URL:
% View data attributes and variables
ncdisp('https://oceanwatch.pifsc.noaa.gov/erddap/griddap/CRW_sst_v3_1_monthly');
Source: https://oceanwatch.pifsc.noaa.gov/erddap/griddap/CRW_sst_v3_1_monthly Format: classic Global Attributes: acknowledgement = 'NOAA Coral Reef Watch program' cdm_data_type = 'Grid' comment = 'This is a product of NOAA Coral Reef Watch Global 5km Satellite Coral Bleaching Heat Stress Monitoring Product Suite Version 3.1, derived from CoralTemp v1.0.' contributor_name = 'NOAA Coral Reef Watch program' contributor_role = 'Collecting source data and deriving products; performing quality control of products; disseminating, storing, and submitting data to archive.' Conventions = 'CF-1.6, ACDD-1.3, COARDS' creator_email = 'coralreefwatch@noaa.gov' creator_institution = 'NOAA/NESDIS/STAR Coral Reef Watch program' creator_name = 'NOAA Coral Reef Watch program' creator_type = 'group' creator_url = 'https://coralreefwatch.noaa.gov' date_created = '2018-03-01T12:00:00Z' date_issued = '2021-03-01T15:05:30Z' date_metadata_modified = '2020-11-18T12:00:00Z' date_modified = '2018-03-01T12:00:00Z' Easternmost_Easting = 359.975 geospatial_bounds = '"POLYGON((-90.0 360.0, 90.0 360.0, 90.0 0.0, -90.0 0.0, -90.0 360.0))"' geospatial_bounds_crs = 'EPSG:4326' geospatial_lat_max = 89.975 geospatial_lat_min = -89.975 geospatial_lat_units = 'degrees_north' geospatial_lon_max = 359.975 geospatial_lon_min = 0.025 geospatial_lon_units = 'degrees_east' grid_mapping_epsg_code = 'EPSG:4326' grid_mapping_inverse_flattening = 298.2572 grid_mapping_name = 'latitude_longitude' grid_mapping_semi_major_axis = 6378137 history = 'Tue Jan 11 11:21:10 2022: ncatted -O -a geospatial_bounds,global,o,c,"POLYGON((-90.0 360.0, 90.0 360.0, 90.0 0.0, -90.0 0.0, -90.0 360.0))" ct5km_sst-mean_v3.1_202112-0-360.nc Tue Jan 11 11:21:10 2022: ncatted -O -a geospatial_lon_max,global,o,f,359.975 ct5km_sst-mean_v3.1_202112-0-360.nc Tue Jan 11 11:21:10 2022: ncatted -O -a geospatial_lon_min,global,o,f,0.025 ct5km_sst-mean_v3.1_202112-0-360.nc Tue Jan 11 11:21:10 2022: ncatted -O -a valid_max,lon,o,f,359.975 ct5km_sst-mean_v3.1_202112-0-360.nc Tue Jan 11 11:21:10 2022: ncatted -O -a valid_min,lon,o,f,0.025 ct5km_sst-mean_v3.1_202112-0-360.nc Tue Jan 11 11:21:06 2022: ncap2 -O -s where(lon<0) lon=lon+360 ct5km_sst-mean_v3.1_202112-0-360.nc ct5km_sst-mean_v3.1_202112-0-360.nc Tue Jan 11 11:21:02 2022: ncks -O --msa_usr_rdr -d lon,0.0,180.0 -d lon,-180.0,0.0 ct5km_sst-mean_v3.1_202112.nc ct5km_sst-mean_v3.1_202112-0-360.nc This is a product data file of the NOAA Coral Reef Watch Global 5km Satellite Coral Bleaching Heat Stress Monitoring Product Suite Version 3.1 (v3.1) in its NetCDF Version 1.0 (v1.0). 2022-04-19T22:05:20Z (local files) 2022-04-19T22:05:20Z https://oceanwatch.pifsc.noaa.gov/erddap/griddap/CRW_sst_v3_1_monthly.das' id = 'Satellite_Global_5km_CoralTemp_SST_Monthly_Mean_Composite' infoUrl = 'https://coralreefwatch.noaa.gov/product/5km/index.php' institution = 'NOAA/NESDIS/STAR Coral Reef Watch program' instrument = 'ATSR-1, ATSR-2, AATSR, AVHRR, AVHRR-2, AVHRR-3, VIIRS, GOES Imager, MTSAT Imager, MTSAT 2 Imager, AHI, ABI, SEVIRI, buoy - moored buoy, buoy - drifting buoy, buoy - TAO buoy, surface seawater intake' instrument_vocabulary = 'NOAA NODC Ocean Archive System Instruments' keywords = '5km, analysed, composite, coral, coraltemp, crw, data, earth, Earth Science > Oceans > Ocean Temperature > Sea Surface Temperature, Earth Science > Oceans > Ocean Temperature > Water Temperature, Earth Science > Spectral/Engineering > Infrared Wavelengths > Thermal Infrared, engineering, environmental, global, information, infrared, mean, month, monthly, national, nesdis, noaa, ocean, oceans, program, reef, satellite, science, sea, sea_surface_temperature, seawater, service, spectral, spectral/engineering, sst, star, surface, temperature, thermal, time, version, watch, water, wavelengths' keywords_vocabulary = 'GCMD Science Keywords' license = 'The data produced by Coral Reef Watch are available for use without restriction, but Coral Reef Watch relies on the ethics and integrity of the user to ensure that the source of the data and products is appropriately cited and credited. When using these data and products, credit and courtesy should be given to NOAA Coral Reef Watch. Please include the appropriate DOI associated with this dataset in the citation. For more information, visit the NOAA Coral Reef Watch website: https://coralreefwatch.noaa.gov. Recommendations for citing and providing credit are provided at https://coralreefwatch.noaa.gov/satellite/docs/recommendations_crw_citation.php. Users are referred to the footer section of the Coral Reef Watch website (https://coralreefwatch.noaa.gov/index.php) for disclaimers, policies, notices pertaining to the use of the data.' metadata_link = 'https://coralreefwatch.noaa.gov/product/5km/index.php' naming_authority = 'gov.noaa.coralreefwatch' NCO = 'netCDF Operators version 4.7.5 (Homepage = http://nco.sf.net, Code = https://github.com/nco/nco)' nco_openmp_thread_number = 1 Northernmost_Northing = 89.975 platform = 'Ships, drifting buoys, moored buoys, TOGA-TAO buoy arrays, GOES-8 satellite, GOES-9 satellite, GOES-10 satellite, GOES-11 satellite, GOES-12 satellite, GOES-13 satellite, GOES-14 satellite, GOES-15 satellite, GOES-16 satellite, MTSAT-1R satellite, MTSAT-2 satellite, Himawari-8 satellite, Meteosat-8 satellite, Meteosat-9 satellite, Meteoset-10 satellite, Meteosat-11 satellite, Suomi NPP, MetOp-A satellite, MetOp-B satellite, NOAA-9 satellite, NOAA-11 satellite, NOAA-12 satellite, NOAA-14 satellite, NOAA-15 satellite, NOAA-16 satellite, NOAA-17 satellite, NOAA-18 satellite, NOAA-19 satellite.' platform_vocabulary = 'NOAA NODC Ocean Archive System Platforms' processing_level = 'Derived from L4 satellite sea surface temperaure analysis' product_version = '3.1' program = 'NOAA Coral Reef Watch program' project = 'NOAA Coral Reef Watch program' publisher_email = 'coralreefwatch@noaa.gov' publisher_institution = 'NOAA/NESDIS/STAR Coral Reef Watch program' publisher_name = 'NOAA Coral Reef Watch program' publisher_type = 'group' publisher_url = 'https://coralreefwatch.noaa.gov' references = 'https://coralreefwatch.noaa.gov/product/5km/index.php and https://coralreefwatch.noaa.gov/satellite/coraltemp.php' source = 'Coral Reef Watch CoralTemp v1.0' sourceUrl = '(local files)' Southernmost_Northing = -89.975 spatial_resolution = '0.05 degree' standard_name_vocabulary = 'CF Standard Name Table v27' summary = 'NOAA Coral Reef Watch Global 5km Satellite Sea Surface Temperature (CoralTemp Version 1.0) Monthly Mean Composite for February 2021. This is a product of NOAA Coral Reef Watch Global 5km Satellite Coral Bleaching Heat Stress Monitoring Product Suite, derived from CoralTemp v1.0.' time_coverage_duration = 'P1M' time_coverage_end = '2021-12-31T12:00:00Z' time_coverage_resolution = 'P1M' time_coverage_start = '1985-01-31T12:00:00Z' title = 'Sea Surface Temperature, Coral Reef Watch, CoralTemp, v3.1 - Monthly, 1985-present' Westernmost_Easting = 0.025 Dimensions: latitude = 3600 longitude = 7200 time = 444 Variables: time Size: 444x1 Dimensions: time Datatype: double Attributes: _CoordinateAxisType = 'Time' actual_range = [476020800 1640952000] axis = 'T' coverage_content_type = 'coordinate' ioos_category = 'Time' long_name = 'reference time of the last day of the composite temporal coverage' standard_name = 'time' time_origin = '01-JAN-1970 00:00:00' units = 'seconds since 1970-01-01T00:00:00Z' latitude Size: 3600x1 Dimensions: latitude Datatype: single Attributes: _CoordinateAxisType = 'Lat' actual_range = [-89.975 89.975] axis = 'Y' comment = 'equirectangular projection and grid centers' coverage_content_type = 'coordinate' ioos_category = 'Location' long_name = 'Latitude' standard_name = 'latitude' units = 'degrees_north' valid_max = 89.975 valid_min = -89.975 longitude Size: 7200x1 Dimensions: longitude Datatype: single Attributes: _CoordinateAxisType = 'Lon' actual_range = [0.025 359.975] axis = 'X' comment = 'equirectangular projection and grid centers' coverage_content_type = 'coordinate' ioos_category = 'Location' long_name = 'Longitude' standard_name = 'longitude' units = 'degrees_east' valid_max = 359.975 valid_min = 0.025 sea_surface_temperature Size: 7200x3600x444 Dimensions: longitude,latitude,time Datatype: single Attributes: _FillValue = NaN colorBarMaximum = 32 colorBarMinimum = 0 coverage_content_type = 'physicalMeasurement' ioos_category = 'Temperature' long_name = 'analysed sea surface temperature' standard_name = 'sea_surface_temperature' units = 'degree_C' valid_max = 50 valid_min = -2
Notice that the full data set has four variables: time, latitude, longitude, and sea_surface_temperature. Rather than downloading all these data, which would be very slow, we'll use this information to access just the data we're interested in. Note: this tutorial uses NetCDF files, but ERDDAP also allows you to download your data of interest as a MATLAB binary file (.mat). A number of other data format options are avalable in the 'File type' drop down shown above.

Time

Notice in the output above that the units for time are 'seconds since 1970-01-01T00:00:00Z'. We can convert this into something that's easier for us to work with. Matlab will convert this type of time to a friendly format using the 'datevec' function. But, Matlab does this assuming that the number you give it represents the number of days since January 0, 0000. We'll get around these differences using the code below to align the two different ways of keeping time by:
  1. Converting 'seconds since...' to 'days since...' by dividing the output time by 86400 (60 seconds in a minute x 60 minutes in an hour x 24 hours in a day).
  2. Adding the time between Jan 0, 0000 and 1970-01-01 using the 'datenum' fuction. 'datenum' does the opposite of 'datevec'. It converts time in [Y M D H M S] to days since January 0, 0000.
And then we'll convert the time to a friendly [Y M D H M S] format using the 'datevec' function.
% Read the time indices, in their native units (seconds since
% 1970-01-01T00:00:00Z')
time_full = ncread('https://oceanwatch.pifsc.noaa.gov/erddap/griddap/CRW_sst_v3_1_monthly', 'time');
% Convert this to [Y M D H M S]
time_full_ymdhms = datevec(time_full/86400 + datenum([1970 1 1 0 0 0]));
Now we have an easy-to-use time index ('time_ymdhms') so we can access just the year we're interested in: 2018.
% Find 2018 (aoi = area of interest)
time_aoi = find(time_full_ymdhms(:,1) == 2018);

Latitude and Longitude

Before we can create a map, we also need to create indices for latitude and longitude. The output above shows us that the units for these are 'degrees_north' and 'degrees_east', repectively.
lat_full = ncread('https://oceanwatch.pifsc.noaa.gov/erddap/griddap/CRW_sst_v3_1_monthly', 'latitude');
lon_full = ncread('https://oceanwatch.pifsc.noaa.gov/erddap/griddap/CRW_sst_v3_1_monthly', 'longitude');
Notice in the output above that sea_surface_temperature has the dimensions of longitude x latitude x time which means it has a size of 7200 x 3600 x 444. In the code below, we'll create indices for our particular area of interest in the central North Pacific.
% Find longitudes from 195 - 210 E
lon_aoi = find(lon_full >= 195 & lon_full <= 210);
% Find latitudes from 17 - 30 N
lat_aoi = find(lat_full >= 17 & lat_full <= 30);

Sea Surface Temperature

Now, we'll access just the small amount of sea_surface_temperature data we're interested in. We'll do this by telling the computer to access the variable 'sea_surface_temperature', beginning at a specific location and time, and spanning our area and time of interest.
% Start coordinates
aoi_start = [lon_aoi(1) lat_aoi(1) time_aoi(1)];
% Coordinates to span
aoi_span = [length(lon_aoi) length(lat_aoi) length(time_aoi)];
sst = ncread('https://oceanwatch.pifsc.noaa.gov/erddap/griddap/CRW_sst_v3_1_monthly', ...
'sea_surface_temperature', aoi_start, aoi_span);
Notice that our workspace now has a sst variable that's sized 300 x 260 x 12. Before we see what these data look like, let's create indices for the area and timespan of data we downloaded and then tidy up our workspace. We'll need the indices when we plot the data.
% Area and time indices, in a format Matlab is expecting
lat = double(lat_full(lat_aoi));
lon = double(lon_full(lon_aoi));
time = time_full(time_aoi);
time_ymdhms = time_full_ymdhms(time_aoi,:); % Note that this variable has 6 columns, unlike the others
% Tidying up
clear *aoi* *full* % Delete every variable with 'aoi' or with 'full' in its name

Working with the extracted data

Let's create a map for a single month, January 2018, the first month in the time span we downloaded.
figure
axesm('mercator', 'MapLatLimit', [17 30], 'MapLonLimit', [195 210], 'MeridianLabel', 'on', ...
'ParallelLabel', 'on', 'MLabelLocation', -164:2:-150, 'PLabelLocation', 18:2:30, ...
'MLabelParallel', 'south');
contourm(lat, lon, sst(:,:,1)', 100, 'fill', 'on');
title(sprintf('Monthly SST %s', datestr(time_ymdhms(1,:), 'mmm-yyyy')));
% Set the color map to jet colors, with 100 levels
colormap(jet(100));
% Add a color bar and label it
c = colorbar;
c.Label.String = 'SST';
% Example of how to mark points on the map
plotm(repelem(26,4), 202:1:205, 'ko', 'MarkerFaceColor', 'k');
% Example of how to add a specific contour, here 20
contourm(lat, lon, sst(:,:,1)', [20 20], 'k');
% Add land and color it grey
geoshow('landareas.shp', 'FaceColor', [0.5 0.5 0.5]);
tightmap % This removes an additional frame around the map that can interfere with the labeling

Plotting a Time Series

Let's pick the following box: 24-26N, 200-206E. We are going to generate a time series of mean SST within that box. To do this, we'll use the same strategy we used above to identify a subsetted area of interest.
% Find longitudes from 200 - 206 E
lon_aoi = find(lon >= 200 & lon <= 206);
% Find latitudes from 24 - 26 N
lat_aoi = find(lat >= 24 & lat <= 26);
% Subset our new area of interest
sst_subset = sst(lon_aoi, lat_aoi,:);
% Average over the area for each of the 12 months
sst_ts(1:12,1) = NaN;
for m = 1:1:12
sst_ts(m,1) = mean(sst_subset(:,:,m), "all", "omitnan");
end
% Plot
figure
plot(time_ymdhms(:,2), sst_ts, 'k-o', 'MarkerFaceColor', 'k');
xlabel('Month');
ylabel('SST (Deg C)');

Creating a map of average SST over a year

We can also create a map of SST averaged of our full time period of interest. Let's go back to using the full area we downloaded, too.
% Average over time, which is the third dimention of our data
sst_yr = mean(sst, 3, "omitnan");
% Plot
figure
axesm('mercator', 'MapLatLimit', [17 30], 'MapLonLimit', [195 210], 'MeridianLabel', 'on', ...
'ParallelLabel', 'on', 'MLabelLocation', -164:2:-150, 'PLabelLocation', 18:2:30, ...
'MLabelParallel', 'south');
contourm(lat, lon, sst_yr', 100, 'fill', 'on');
title([sprintf('Mean SST %s', datestr(time_ymdhms(1,:),'mmm-yyyy')) ...
sprintf(' - %s', datestr(time_ymdhms(12,:),'mmm-yyyy'))]);
% Set the color map to jet colors, with 100 levels
colormap(jet(100));
% Add a color bar and label it
c = colorbar;
c.Label.String = 'SST';
% Add land and color it grey
geoshow('landareas.shp', 'FaceColor', [0.5 0.5 0.5]);
tightmap