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 daily SST data:
https://coastwatch.glerl.noaa.gov/erddap/griddap/GLSEA_GCS.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://coastwatch.glerl.noaa.gov/erddap/griddap/GLSEA_GCS.nc?sst%5B(2021-07-21T12:00:00Z):1:(2021-07-28T12:00:00Z)%5D%5B(38.8749871947229):1:(50.6059751976437)%5D%5B(-92.4199507342304):1:(-75.8816402880577)%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://coastwatch.glerl.noaa.gov/erddap/griddap/GLSEA_GCS');
Source: https://coastwatch.glerl.noaa.gov/erddap/griddap/GLSEA_GCS Format: classic Global Attributes: cdm_data_type = 'Grid' Conventions = 'CF-1.6, COARDS, ACDD-1.3' Easternmost_Easting = -75.8816 GDAL = 'GDAL 2.1.3, released 2017/20/01' geospatial_lat_max = 50.606 geospatial_lat_min = 38.875 geospatial_lat_resolution = 0.014016 geospatial_lat_units = 'degrees_north' geospatial_lon_max = -75.8816 geospatial_lon_min = -92.42 geospatial_lon_resolution = 0.014016 geospatial_lon_units = 'degrees_east' history = 'Sea Surface Temperature (SST) from Great Lakes Surface Environmental Analysis (GLSEA) asc format to nc fromat 2022-04-15T18:56:18Z (local files) 2022-04-15T18:56:18Z http://coastwatch.glerl.noaa.gov/erddap/griddap/GLSEA_GCS.das' infoUrl = 'https://coastwatch.glerl.noaa.gov/glsea/glsea.html' institution = 'CoastWatch Great Lakes Node' keywords = 'analysis, data, earth, Earth Science > Oceans > Ocean Temperature > Water Temperature, environmental, glsea, great, great lakes, lakes, ocean, oceans, science, sea, sea_water_temperature, seawater, sst, surface, temperature, time, water' keywords_vocabulary = 'GCMD Science Keywords' license = 'The data may be used and redistributed for free but is not intended for legal use, since it may contain inaccuracies. Neither the data Contributor, ERD, NOAA, nor the United States Government, nor any of their employees or contractors, makes any warranty, express or implied, including warranties of merchantability and fitness for a particular purpose, or assumes any legal liability for the accuracy, completeness, or usefulness, of this information.' NCO = 'netCDF Operators version 4.7.5 (Homepage = http://nco.sf.net, Code = https://github.com/nco/nco)' Northernmost_Northing = 50.606 source = 'Sea Surface Temperature (SST) from Great Lakes Surface Environmental Analysis (GLSEA)' sourceUrl = '(local files)' Southernmost_Northing = 38.875 standard_name_vocabulary = 'CF Standard Name Table v29' summary = 'Sea Surface Temperature (SST) from Great Lakes Surface Environmental Analysis (GLSEA)' time_coverage_end = '2022-04-14T12:00:00Z' time_coverage_start = '1995-01-01T12:00:00Z' title = 'Sea Surface Temperature (SST) from Great Lakes Surface Environmental Analysis (GLSEA), Geodetic coordinate system (LAT, LON), 1995-present' Westernmost_Easting = -92.42 Dimensions: latitude = 838 longitude = 1181 time = 9962 Variables: time Size: 9962x1 Dimensions: time Datatype: double Attributes: _CoordinateAxisType = 'Time' actual_range = [788961600 1649937600] axis = 'T' calendar = 'Gregorian' ioos_category = 'Time' long_name = 'Easten Time' standard_name = 'time' time_origin = '01-JAN-1970 00:00:00' units = 'seconds since 1970-01-01T00:00:00Z' latitude Size: 838x1 Dimensions: latitude Datatype: double Attributes: _CoordinateAxisType = 'Lat' actual_range = [38.875 50.606] axis = 'Y' ioos_category = 'Location' long_name = 'Latitude' standard_name = 'latitude' units = 'degrees_north' longitude Size: 1181x1 Dimensions: longitude Datatype: double Attributes: _CoordinateAxisType = 'Lon' actual_range = [-92.42 -75.8816] axis = 'X' ioos_category = 'Location' long_name = 'Longitude' standard_name = 'longitude' units = 'degrees_east' sst Size: 1181x838x9962 Dimensions: longitude,latitude,time Datatype: single Attributes: _FillValue = -99999 colorBarMaximum = 32 colorBarMinimum = 0 grid_mapping = 'crs' ioos_category = 'Temperature' long_name = 'Temperature' standard_name = 'sea_water_temperature' units = 'degree_C'
Notice that the full data set has four variables: time, latitude, longitude, and sst. 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://coastwatch.glerl.noaa.gov/erddap/griddap/GLSEA_GCS', '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 week we're interested in: 21 - 28 July 2021.
% Find 2018 (aoi = area of interest)
time_aoi = find(time_full_ymdhms(:,1) == 2021 & time_full_ymdhms(:,2) == 7 ...
& time_full_ymdhms(:,3) >= 21 & time_full_ymdhms(:,3) <= 28);

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 = ncread('https://coastwatch.glerl.noaa.gov/erddap/griddap/GLSEA_GCS', 'latitude');
lon = ncread('https://coastwatch.glerl.noaa.gov/erddap/griddap/GLSEA_GCS', 'longitude');
Notice in the output above that sst has the dimensions of longitude x latitude x time which means it has a size of 1181 x 838 x 9961.

Sea Surface Temperature

Now, we'll access just the small amount of sst data we're interested in. We'll do this by telling the computer to access the variable 'sst', beginning at a specific time and spanning our time of interest.
% Start coordinates
aoi_start = [1 1 time_aoi(1)];
% Coordinates to span
aoi_span = [length(lon) length(lat) length(time_aoi)];
sst = ncread('https://coastwatch.glerl.noaa.gov/erddap/griddap/GLSEA_GCS', 'sst', aoi_start, aoi_span);
Notice that our workspace now has a sst variable that's sized 1181 x 838 x 8. Before we see what these data look like, let's create indices for the 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
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* % Delete every variable with 'aoi' in its name

Working with the extracted data

Let's create a map for a single day, 21 July 2021, the first day in the time span we downloaded.
figure
axesm('mercator', 'MapLatLimit', [38.8 50.7], 'MapLonLimit', [-92.4 -75.8], 'MeridianLabel', 'on', ...
'ParallelLabel', 'on', 'MLabelLocation', -90:5:-80, 'PLabelLocation', 40:5:50, ...
'MLabelParallel', 'south');
contourm(lat, lon, sst(:,:,1)', 10:0.5:25.5, 'fill', 'on');
title(sprintf('Daily SST %s', datestr(time_ymdhms(1,:), 'dd-mmm-yyyy')));
% Set the color map to jet colors, with 100 levels
colormap(jet(32));
% Add a color bar and label it
c = colorbar;
c.Label.String = 'SST';
% Example of how to mark points on the map
plotm(repelem(45,5), -83:0.5:-81, 'ko', 'MarkerFaceColor', 'k');
% Example of how to add a specific contour, here 15
contourm(lat, lon, sst(:,:,1)', [15 15], 'k');
% Add land and color it grey
geoshow('landareas.shp', 'FaceColor', [0.5 0.5 0.5]);
% Example of how to add a marker to tutorial author's home town, just for
% fun
plotm(45.8527, -87.0218, 'p');
tightmap % This removes an additional frame around the map that can interfere with the labeling

Plotting a Time Series

Let's pick a box encompassing Lake Superior: north of 46N and west of 84W. 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 for time to identify a subsetted area of interest.
% Find longitudes west of 84 W
lon_aoi = find(lon <= -84);
% Find latitudes north of 46 N
lat_aoi = find(lat >= 46);
% Subset our new area of interest
sst_subset = sst(lon_aoi, lat_aoi,:);
% Average over the area for each of the 8 days
sst_ts(1:8,1) = NaN;
for d = 1:1:8
sst_ts(d,1) = mean(sst_subset(:,:,d), "all", "omitnan");
end
% Plot
figure
plot(time_ymdhms(:,3), sst_ts, 'k-o', 'MarkerFaceColor', 'k');
xlabel('Month');
ylabel('SST (Deg C)');
set(gca, 'XTick', time_ymdhms(:,3));
set(gca, 'XTickLabel', datestr(time_ymdhms(:,:), 'mm/dd'))

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_wk = mean(sst, 3, "omitnan");
figure
axesm('mercator', 'MapLatLimit', [38.8 50.7], 'MapLonLimit', [-92.4 -75.8], 'MeridianLabel', 'on', ...
'ParallelLabel', 'on', 'MLabelLocation', -90:5:-80, 'PLabelLocation', 40:5:50, ...
'MLabelParallel', 'south');
contourm(lat, lon, sst_wk', 10:0.5:25.5, 'fill', 'on');
title([sprintf('Mean SST %s', datestr(time_ymdhms(1,:),'dd-mmm-yyyy')) ...
sprintf(' - %s', datestr(time_ymdhms(8,:),'dd-mmm-yyyy'))]);
% Set the color map to jet colors, with 100 levels
colormap(jet(32));
% 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