% Project: Assessing storm energy reduction by the vegetated salt marsh platform in Newbury, MA, pressure timeseries
% Sample code for accessing the pressure timeseries data (DOI:10.5281/zenodo.15635830) 
% Created by M. Foster-Martinez (6/16/2025)

%User input
stationID = 'ETF'; %Nine different sites, each with a three letter code
path = ['F:\A_drive_Oct2021\Newbury_original_data' filesep stationID];

%Retrieve and read netCDF
fname_dum = dir([path, filesep,'NB18',stationID,'rbr.cdf']);
nc_fname = fname_dum.name;
cd(path)
info_nc = ncinfo(nc_fname);
warning('off','all');

%Attribute information to be used in calculations
ht_m = ncreadatt(nc_fname,'/','initial_instrument_height'); %meters
deploy_t = ncreadatt(nc_fname,'/','Deployment_date');
recovery_t = ncreadatt(nc_fname,'/','Recovery_date');
BurstSampleCount = ncreadatt(nc_fname,'/','samples_per_burst'); %count
LoggingSampling =round(1/ ncreadatt(nc_fname,'/','sample_interval')); %Hz

% time and pressure data 
tim1 = double(ncread(nc_fname,'time'));
tim2 = double(ncread(nc_fname,'time2'));
jd_big = tim1+tim2/86400000;
if strcmp(stationID,'ETF')
    jd_big = jd_big - 0.25; %EFT was on UTC time instead of eastern. Eastern = UTC - 6 hr. 6 hr is 0.25 of one day. so subtract a quarter of a day to get this to eastern time. 
end
jd = jd_big(1,~isnan(jd_big(1,:)));
indnan = ~isnan(jd);
dat_db = double(ncread(nc_fname,'Pressure')); %dbar
if strcmp(stationID,'STF')
    dat_db = dat_db(:,1:2332); %cut filler values in this matrix 
end
deploy_jd = date2jd(deploy_t);
recovery_jd = date2jd(recovery_t);
inWater = find(jd>deploy_jd & jd<recovery_jd);