3. Extract data within a shapefile using ERDDAP

This tutorial uses the Xtractomatic package to download data from within the boundaries of the Papahanaumokuakea Marine National Monument (PMNM).

This tutorial will teach you how to extract and display SST values for a particular time period or average SST over the whole time-series available within a shapefile.

The shapefile for the PMNM boundaries can be downloaded here: http://sanctuaries.noaa.gov/library/imast_gis.html. Save the file https://sanctuaries.noaa.gov/library/imast/pmnm_py.zip on your computer and extract it. It would be easiest for the purposes of this tutorial if you place the pmnm_py folder in the directory you're working in. But, if you want to keep it elsewhere, you'll just need to edit the path name accordingly.

Data extraction

We'll do this in two steps. First, we'll download the satellite data for a footprint that's slightly larger than our area of interest. Then, we'll clip these data to only those which fall within the Monument.

Download satellite data of interest

% As always, it's helpful it take a look at what we're working with: 
ncdisp('https://oceanwatch.pifsc.noaa.gov/erddap/griddap/goes-poes-monthly-ghrsst-RAN');
Source:
           https://oceanwatch.pifsc.noaa.gov/erddap/griddap/goes-poes-monthly-ghrsst-RAN
Format:
           classic
Global Attributes:
           acknowledgement           = 'NOAA/NESDIS'
           cdm_data_type             = 'Grid'
           comment                   = 'The Geo-Polar Blended Sea Surface Temperature (SST) Analysis combines multi-satellite retrievals of sea surface temperature into a single analysis of SST. This analysis includes only nighttime data.'
           Conventions               = 'CF-1.6, Unidata Observation Dataset v1.0, COARDS, ACDD-1.3'
           creator_email             = 'john.sapper@noaa.gov'
           creator_name              = 'Office of Satellite Products and Operations'
           creator_type              = 'group'
           creator_url               = 'www.osdpd.nesdis.noaa.gov'
           date_created              = '2018-12-26T10:43:22Z'
           Easternmost_Easting       = 359.975
           gds_version_id            = '2.0'
           geospatial_lat_max        = 89.975
           geospatial_lat_min        = -89.975
           geospatial_lat_resolution = 0.05
           geospatial_lat_units      = 'degrees_north'
           geospatial_lon_max        = 359.975
           geospatial_lon_min        = 0.025
           geospatial_lon_resolution = 0.05
           geospatial_lon_units      = 'degrees_east'
           history                   = 'Fri Apr 15 05:39:25 2022: ncea -v analysed_sst,sea_ice_fraction /mnt/r01/data/goes-poes_ghrsst/daily/20220301000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220302000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220303000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220304000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220305000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220306000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220307000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220308000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220309000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220310000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220311000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220312000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220313000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220314000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220315000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220316000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220317000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220318000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220319000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220320000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220321000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220322000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220323000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220324000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220325000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220326000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220327000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220328000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220329000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220330000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc /mnt/r01/data/goes-poes_ghrsst/daily/20220331000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc 202203-gp-monthly.nc
                                       Wed Mar  2 02:05:24 2022: ncatted -O -a easternmost_longitude,global,o,f,360.0 20220301000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc
                                       Wed Mar  2 02:05:24 2022: ncatted -O -a westernmost_longitude,global,o,f,0.0 20220301000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc
                                       Wed Mar  2 02:05:24 2022: ncatted -O -a valid_max,lon,o,f,360.0 20220301000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc
                                       Wed Mar  2 02:05:24 2022: ncatted -O -a valid_min,lon,o,f,0.0 20220301000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc
                                       Wed Mar  2 02:05:16 2022: ncap2 -O -s where(lon<0) lon=lon+360 20220301000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc 20220301000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc
                                       Wed Mar  2 02:05:08 2022: ncks -O --msa_usr_rdr -d lon,0.0,180.0 -d lon,-180.0,0.0 20220301000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0.nc 20220301000000-OSPO-L4_GHRSST-SSTfnd-Geo_Polar_Blended_Night-GLOB-v02.0-fv01.0-0-360.nc
                                       NESDIS geo-SST L1 to L2 processor, NESDIS Advanced Clear-Sky Processor for Oceans (ACSPO), NESDIS Geo-Polar 1/20th degree Blended SST Analysis
                                       2022-04-15T23:55:19Z (local files)
                                       2022-04-15T23:55:19Z https://oceanwatch.pifsc.noaa.gov/erddap/griddap/goes-poes-monthly-ghrsst-RAN.das'
           id                        = 'Geo_Polar_Blended_Night-OSPO-L4-GLOB-v1.0'
           infoUrl                   = 'https://podaac.jpl.nasa.gov/dataset/Geo_Polar_Blended_Night-OSPO-L4-GLOB-v1.0'
           institution               = 'Office of Satellite Products and Operations'
           keywords                  = 'analysed, analysed_sst, area, blended, cryosphere, data, distribution, earth, Earth Science > Cryosphere > Sea Ice > Ice Extent, Earth Science > Oceans > Ocean Temperature > Sea Surface Temperature, Earth Science > Oceans > Sea Ice > Ice Extent, extent, foundation, fraction, global, ice, ice distribution, infrared, input, night, ocean, oceans, office, only, operations, over, products, satellite, science, sea, sea_ice_area_fraction, sea_ice_fraction, sea_surface_foundation_temperature, sst, surface, temperature, time, using'
           keywords_vocabulary       = 'GCMD Science Keywords'
           license                   = 'GHRSST protocol describes data use as free and open'
           metadata_link             = 'https://podaac.jpl.nasa.gov/ws/metadata/dataset?format=iso&shortName=Geo_Polar_Blended_Night-OSPO-L4-GLOB-v1.0'
           naming_authority          = 'org.ghrsst'
           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                  = 'Suomi NPP, MetOp-B, GOESE (GOES-16), GOESW (GOES-15), Meteosat-11, Meteosat-08, Himawari-8'
           processing_level          = 'L4'
           product_version           = '1.0'
           project                   = 'Group for High Resolution Sea Surface Temperature'
           publisher_email           = 'ghrsst-po@nceo.ac.uk'
           publisher_name            = 'The GHRSST Project Office'
           publisher_type            = 'group'
           publisher_url             = 'https://www.ghrsst.org'
           references                = 'Fieguth,P.W. et al. "Mapping Mediterranean altimeter data with a multiresolution optimal interpolation algorithm", J. Atmos. Ocean Tech, 15 (2): 535-546, 1998.     Fieguth, P. Multiply-Rooted Multiscale Models for Large-Scale Estimation, IEEE Image Processing, 10(11), 1676-1686, 2001.     Khellah, F., P.W. Fieguth, M.J. Murray and M.R. Allen, "Statistical Processing of Large Image Sequences", IEEE Transactions on Geoscience and Remote Sensing, 12 (1), 80-93, 2005.'
           sensor                    = 'VIIRS, AVHRR, ABI, Imager, SEVIRI, SEVIRI, AHI'
           source                    = 'OSPO-ACSPO_VIIRS, OSPO-ACSPO_METOPB_FRAC, OSPO-GOES16_SST_L3, OSPO-GOES15_SST_L3, OSPO-METEOSAT11_SST_L3, OSPO-HIMA8_SST_L3'
           sourceUrl                 = '(local files)'
           Southernmost_Northing     = -89.975
           spatial_resolution        = '0.05 degree'
           standard_name_vocabulary  = 'CF Standard Name Table v29'
           summary                   = 'Monthly composite of analysed blended sea surface temperature over the global ocean using night only input data. An SST estimation scheme which combines multi-satellite retrievals of sea surface temperature datasets available from polar orbiters, geostationary InfraRed (IR) and microwave sensors into a single global analysis. This global SST analysis provides a daily gap free map of the foundation sea surface temperature at 0.05� spatial resolution. The monthly composites were generated from daily files by OceanWatch Central Pacific'
           time_coverage_end         = '2022-03-01T12:00:00Z'
           time_coverage_start       = '2002-09-05T12:00:00Z'
           title                     = 'Sea Surface Temperature, NOAA geopolar blended - Monthly, 2002-Present (2017 Reanalysis)'
           Westernmost_Easting       = 0.025
Dimensions:
           latitude  = 3600
           longitude = 7200
           time      = 227
Variables:
    time            
           Size:       227x1
           Dimensions: time
           Datatype:   double
           Attributes:
                       _CoordinateAxisType = 'Time'
                       actual_range        = [1031227200  1646136000]
                       axis                = 'T'
                       calendar            = 'Gregorian'
                       comment             = 'Nominal time of Level 4 analysis'
                       ioos_category       = 'Time'
                       long_name           = 'reference time of sst field'
                       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'
                       ioos_category       = 'Location'
                       long_name           = 'Latitude'
                       standard_name       = 'latitude'
                       units               = 'degrees_north'
                       valid_max           = 90
                       valid_min           = -90
    longitude       
           Size:       7200x1
           Dimensions: longitude
           Datatype:   single
           Attributes:
                       _CoordinateAxisType = 'Lon'
                       actual_range        = [0.025       359.975]
                       axis                = 'X'
                       comment             = 'equirectangular projection'
                       ioos_category       = 'Location'
                       long_name           = 'Longitude'
                       standard_name       = 'longitude'
                       units               = 'degrees_east'
                       valid_max           = 360
                       valid_min           = 0
    analysed_sst    
           Size:       7200x3600x227
           Dimensions: longitude,latitude,time
           Datatype:   single
           Attributes:
                       _FillValue      = NaN
                       colorBarMaximum = 305
                       colorBarMinimum = 273
                       comment         = 'nighttime analysed SST for each ocean grid point'
                       ioos_category   = 'Temperature'
                       long_name       = 'analysed sea surface temperature'
                       references      = 'Fieguth,P.W. et al. "Mapping Mediterranean altimeter data with a multiresolution optimal interpolation algorithm", J. Atmos. Ocean Tech, 15 (2): 535-546, 1998.     Fieguth, P. Multiply-Rooted Multiscale Models for Large-Scale Estimation, IEEE Image Processing, 10(11), 1676-1686, 2001.     Khellah, F., P.W. Fieguth, M.J. Murray and M.R. Allen, "Statistical Processing of Large Image Sequences", IEEE Transactions on Geoscience and Remote Sensing, 12 (1), 80-93, 2005.'
                       source          = 'OSPO-ACSPO_VIIRS, OSPO-ACSPO_METOPB_FRAC, OSPO-GOES16_SST_L3,OSPO-GOES15_SST_L3, OSPO-METEOSAT11_SST_L3, OSPO-HIMA8_SST_L3'
                       standard_name   = 'sea_surface_foundation_temperature'
                       units           = 'degree_K'
                       valid_max       = 4000
                       valid_min       = -200
    sea_ice_fraction
           Size:       7200x3600x227
           Dimensions: longitude,latitude,time
           Datatype:   single
           Attributes:
                       _FillValue      = NaN
                       colorBarMaximum = 1
                       colorBarMinimum = 0
                       comment         = ' Percentage of ice'
                       ioos_category   = 'Ice Distribution'
                       long_name       = 'sea ice fraction'
                       source          = 'NCEP 1/12th degree ice mask'
                       standard_name   = 'sea_ice_area_fraction'
                       units           = '1'
                       valid_max       = 100
                       valid_min       = 0

We can see in the information above that longitude is in the 0 - 360 degree format. And we know from reading the ReadMe file for the monument boundaries that the shapefile longitude is in the -180 - 180 degree format. This means we need convert the shapefile longitudes to a 0 - 360 format for indexing puposes (the native format will plot just fine).

% Access the Lat and Lon coordinates
PMNM = shaperead('pmnm_py/PMNM_py_files/PMNM_py.shp','UseGeoCoords',true);
% Convert negative longitudes to positive values
neg_lon = find(PMNM.Lon < 0);
PMNM.Lon(neg_lon) = PMNM.Lon(neg_lon) + 360;
clear neg_lon
% Find the minimum and maximum lat and lon values, for using to download
% SST data
[Lon_min, Lon_max] = bounds(PMNM.Lon);
[Lat_min, Lat_max] = bounds(PMNM.Lat);
% Now we can follow the steps we used in the earlier tutorials to download
% the SST data
time_full = ncread('https://oceanwatch.pifsc.noaa.gov/erddap/griddap/goes-poes-monthly-ghrsst-RAN', ...
   'time');
% Convert this to [Y M D H M S]
time_full_ymdhms = datevec(time_full/86400 + datenum([1970 1 1 0 0 0])); 
% In this particular example, we're interested in March - November 2015
time_aoi = find(time_full_ymdhms(:,1) == 2015 & time_full_ymdhms(:,2) >= 3 & ...
    time_full_ymdhms(:,2) <= 11);
lat_full = ncread('https://oceanwatch.pifsc.noaa.gov/erddap/griddap/goes-poes-monthly-ghrsst-RAN', ...
    'latitude');
lon_full = ncread('https://oceanwatch.pifsc.noaa.gov/erddap/griddap/goes-poes-monthly-ghrsst-RAN', ...
    'longitude');
% Find longitudes that span our area of interest
% Rounding the minimum down and the maximum up
lon_aoi = find(lon_full >= floor(Lon_min) & lon_full <= ceil(Lon_max));
% Find latitudes that span our area of interest
% Rounding the minimum down and the maximum up
lat_aoi = find(lat_full >= floor(Lat_min) & lat_full <= ceil(Lat_max));
% 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)];
% Download the data of interest
sst = ncread('https://oceanwatch.pifsc.noaa.gov/erddap/griddap/goes-poes-monthly-ghrsst-RAN', ...
    'analysed_sst', aoi_start, aoi_span);
% 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* *min *max

Clip the data to fit the Monument boundaries

Now we have everything we need to clip the SST data to just those within the boundaries of the Monument. The code below borrows from the no-longer-supported xtractoMatlab: https://github.com/rmendels/xtractoMatlab. You can learn about other avenues to use this tool at: https://coastwatch.pfeg.noaa.gov/xtracto/.

% The general premise here is that we're going to create a mesh grid of our
% full domain (what we downloaded), identify all the points that are within
% or on the polygon that defines the Monument, and set to NaN those that
% are not.
% Make the meshgrid and the mask
[XLON, XLAT] = meshgrid(lon, lat);
[IN ON] = inpolygon(XLON, XLAT, PMNM.Lon, PMNM.Lat);
mask2D = IN | ON;
% Replicate it over the number of time steps, which is the third dimension
% of our sst variable
mask3D = permute(repmat(mask2D,[1 1 size(sst,3)]),[2 1 3]);
% Set to NaN all points not within the polygon
sst(~mask3D) = NaN;

Plotting the data

The extracted data contains several time steps (months) of sst data in the monument boundaries. Let's make a plot of the second time step. Also, we can see above where we used 'ncdisp', that the units for these data are Kelvin. Let's changes this to degrees Celsius when we make our map.

figure

axesm('mercator', 'MapLatLimit', [18 32], 'MapLonLimit', [177 207], 'MeridianLabel', 'on', ...

'ParallelLabel', 'on', 'MLabelLocation', -180:5:-155, 'PLabelLocation', 20:5:30, ...

'MLabelParallel', 'south'); % This is the basemap

% Plot SST for the second time step, in degrees C, with 50 contour levels

contourm(lat, lon, sst(:,:,2)'-273.15, 50, 'fill', 'on');

% Title the map

title(sprintf('SST %s', datestr(time_ymdhms(2,:), 'mmm-yyyy')));

% Set the color map to jet colors, with 50 levels

colormap(jet(50));

% 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

On your own!

Plot the average SST for the period we downloaded. Here's a hint to get you started:

sst_avg = mean(sst, 3, "omitnan");