%% Get Blue Marble map of the earth
nasa = wmsfind('nasa', 'SearchField', 'serverurl');

layer = nasa.refine('bluemarbleng', ...
'SearchField', 'layername', 'MatchType', 'exact');

layer = wmsupdate(layer);

%% Set the latitude and longitude limits of our map
latlim = [-90 90];
lonlim = [-90 90];

% Use wmsread to read in the map as an image with reference coordinates
[blueMarble, Rmarble] = wmsread(layer,'Latlim',latlim,'Lonlim',lonlim);

% Show the map
geoshow(blueMarble,Rmarble)

%% Let's look at the NOAA Central Pacific OpenWatch Thredds server
baseurl = 'http://oceanwatch.pifsc.noaa.gov/thredds';
web(baseurl,'-browser');

%% Get Sea Surface Temperature (2-day) from OceanWatch
wmsurl = 'http://oceanwatch.pifsc.noaa.gov/thredds/wms/goes-poes/2day?service=WMS&version=1.3.0&request=GetCapabilities';
info = wmsinfo(wmsurl);
[goes,Rgoes] = wmsread(info.Layer,...
    'ImageHeight',size(blueMarble,1),...
    'ImageWidth',size(blueMarble,2),...
    'Latlim',latlim,'Lonlim',lonlim);
geoshow(goes,Rgoes);

%% Composite the GOES sea surface temperature on top of the map
threshold = 180;
index = any(goes < threshold, 3);
index = cat(3, index, index, index);

combined = blueMarble;
combined(index) = goes(index);

%% Show the compositied image
geoshow(combined,Rmarble);

%% Use the NetCDF to get detailed info about the surface temperature
ncurl = 'http://oceanwatch.pifsc.noaa.gov/thredds/dodsC/goes-poes/2day';
ncdisp(ncurl);

%% Read 'lat','lon', and 'time' variables from this dataset
lat = ncread(ncurl,'lat');
lon = ncread(ncurl,'lon');
time = ncread(ncurl,'time');

%% Get theSurface temperature history of a point just off the coast of NC!
diamondShoals = [35 360-75];
latInd = find(lat == diamondShoals(1));
lonInd = find(lon == diamondShoals(2));

%% Read the surface temperature off the coast of Cape Hatteras
surfaceTemp = ncread(ncurl, 'sst',[lonInd latInd 1],[1 1 319]);
surfaceTemp = reshape(surfaceTemp,1,size(surfaceTemp,3));

%% Recompute the Time relative to Jan-1-0000
timeStamps = time(1:end-1) + datenum('1800-1-1 00:00:00');

%% Show the 2-day average surface temperature for the past year
figure;
plot(timeStamps, surfaceTemp);

% Format and Label the Plot
title('Sea Surface Temperature at Diamond Shoals','FontSize',20);
xlabel('Date','FontSize',16);
ylabel('Temperature (Degrees Celsius)','FontSize',16);
datetick('x','mm/dd/yy');