% $Id$
% $Date$
% $Author$
% $Revision$
% $HeadURL$

% original data
url0 = 'http://dtvirt5.deltares.nl:8080/thredds/dodsC/opendap/rijkswaterstaat/jarkus/profiles/transect_r2013.nc';
% filled data
url = 'http://dtvirt5.deltares.nl:8080/thredds/dodsC/opendap/rijkswaterstaat/jarkus/profiles/transect_r2013_filled.nc';
% target file name
ncfile = 'duros_active_profile_volume.nc';

svninfo = svn_info([mfilename('fullpath') '.m']);
[~, svnversion] = system(['svnversion ' svninfo.path]);

info = ncinfo(url);
tr = jarkus_transects(...
    'url', url,...
    'output', {'id', 'time'});

% time info
years = cellfun(@(dv) str2double(datestr(dv, 'yyyy')), num2cell(nc_cf_time(url, 'time')));
time_coverage_duration = max(years) - min(years);
time_coverage_resolution = min(diff(years));

n = length(tr.id);
idx = true(1,n);
for i = 1:n
    multiWaitbar(...
        'Getting boundary conditions', i/n, ...
        'label', sprintf('transect %i', tr.id(i)));
    bc = get_hr_boundaryconditions(...
        'id', tr.id(i),...
        'datadir', '../data');
    %disp(bc)
    if isempty(bc)
        idx(i) = false;
    else
        idx(i) = all(cellfun(@isfinite, struct2cell(bc)));
    end
end

%%
asid = load('heijer2012_transects.txt');
idx2 = ismember(tr.id, asid) | tr.id >= 1e7;

id = tr.id(idx&idx2);
%id = id(1:1000:end);

n = length(id);
nt = length(tr.time);


dimstruct        = nccreateDimstruct('Name','alongshore','Length',n);
dimstruct(end+1) = nccreateDimstruct('Name','time','Length', nt);
dimstruct(end+1) = nccreateDimstruct('Name','bounds','Length', 2);

varstruct        = nccreateVarstruct('Name', 'time',...
    'Dimensions', {'time'},...
    'Attributes', {'standard_name', 'time',...
		'axis', 'T',...
		'units', 'days since 1970-01-01'});
varstruct(end+1) = nccreateVarstruct('Name', 'id',...
    'Dimensions', {'alongshore'},...
    'Attributes', {'long_name', 'identifier', 'comment', 'sum of area code (x1e6) and alongshore coordinate'});
varstruct(end+1) = nccreateVarstruct('Name','volume',...
    'Dimensions',{'alongshore', 'time'},...
    'Attributes',{'comment', 'volume between NAP-55m and the active part of profile; active part differs per year','units','m2'});
varstruct(end+1) = nccreateVarstruct('Name','mean_volume',...
    'Dimensions',{'alongshore'},...
    'Attributes',{'comment', 'mean volume between NAP-55m and the active part of profile averaged over all years; active part differs per year','units','m2'});
varstruct(end+1) = nccreateVarstruct('Name','volumeC',...
    'Dimensions',{'alongshore', 'time'},...
    'Attributes',{'comment', 'volume between NAP-55m and the active part of profile; active part based on outer boundaries of whole time series','units','m2'});
varstruct(end+1) = nccreateVarstruct('Name','mean_volumeC',...
    'Dimensions',{'alongshore'},...
    'Attributes',{'comment', 'mean volume between NAP-55m and the active part of profile averaged over all years; active part based on outer boundaries of whole time series','units','m2'});
varstruct(end+1) = nccreateVarstruct('Name','bounds',...
    'Dimensions',{'alongshore', 'time', 'bounds'},...
    'Attributes',{'comment', 'bounds of active profile w.r.t. RSP; active part differs per year','units','m'});
varstruct(end+1) = nccreateVarstruct('Name','boundsC',...
    'Dimensions',{'alongshore', 'bounds'},...
    'Attributes',{'comment', 'bounds of active profile w.r.t. RSP; active part based on outer boundaries of whole time series','units','m'});
varstruct(end+1) = nccreateVarstruct('Name','boundsC_lon',...
    'Dimensions',{'alongshore', 'bounds'},...
    'Attributes',{'long_name', 'longitude', 'comment', 'longitude of bounds of active in WGS84; active part based on outer boundaries of whole time series','units','degrees east'});
varstruct(end+1) = nccreateVarstruct('Name','boundsC_lat',...
    'Dimensions',{'alongshore', 'bounds'},...
    'Attributes',{'long_name', 'latitude', 'comment', 'latitude of bounds of active in WGS84; active part based on outer boundaries of whole time series','units','degrees north'});
varstruct(end+1) = nccreateVarstruct('Name','contourdepth',...
    'Dimensions',{'alongshore', 'time'},...
    'Attributes',{'comment', 'depth contour at seaward boundary of active part; active part differs per year','units','m'});
varstruct(end+1) = nccreateVarstruct('Name','contourdepthC',...
    'Dimensions',{'alongshore', 'time'},...
    'Attributes',{'comment', 'depth contour at seaward boundary of active part; active part based on outer boundaries of whole time series','units','m'});

datefmt = 'yyyy-mm-ddTHH:MMZ';
jrkid = info.Attributes(ismember({info.Attributes.Name}, 'id')).Value;
% put these in a schema
schema = nccreateSchema(dimstruct, varstruct,...
    'Attributes', {
    'naming_authority','deltares.nl',...
    'id', sprintf('duros_active_profile_volume_%s', datestr(nowutc, 'yyyymmdd')),...
    'Metadata_Conventions', 'Unidata Dataset Discovery v1.0',...
    'title', 'volumes and boundaries of DUROS+ active profile part',...
    'summary', 'volumes below the active sections of the cross-shore transects due to a normative dune erosion event as calculated by DUROS+',...
    'keywords', 'dune erosion, volume, duros+',...
    'history', sprintf('%s: volumes and boundaries calculated based on JarKus transects with id %s using %s (revision %s)', datestr(nowutc, datefmt), jrkid, svninfo.url, strtrim(svnversion)),...
    'comment', 'the NAP-55m level is chosen as lower boundary for the profile volumes because this is just below the lowest measured point in the whole dataset',...
    ...
    'geospatial_lat_min', [],...
    'geospatial_lat_max', [],...
    'geospatial_lon_min', [],...
    'geospatial_lon_max', [],...
    'time_coverage_start', min(years),...
    'time_coverage_end', max(years),...
    'geospatial_vertical_min', [],...
    'geospatial_vertical_max', [],...
    ...
    'geospatial_lon_units', 'decimal degrees east',...
    'geospatial_lon_resolution', [],...
    'geospatial_lat_units', 'decimal degrees north',...
    'geospatial_lat_resolution', [],...
    'geospatial_vertical_units', 'm',...
    'geospatial_vertical_resolution', [],...
    'geospatial_vertical_positive', 'up',...
    'time_coverage_units', 'year',...
    'time_coverage_duration', time_coverage_duration,...
    'time_coverage_resolution', time_coverage_resolution,...
    ...
    'creator_name', 'Kees den Heijer',...
    'creator_url', 'http://www.deltares.nl',...
    'creator_email', 'kees.denheijer@deltares.nl',...
    'institution', 'Deltares',...
    'date_created', datestr(nowutc, datefmt),...
    'date_modified', datestr(nowutc, datefmt),...
    'date_issued', datestr(nowutc, datefmt),...
    'project', 'WTI dunes temporal effects',...
...    'acknowledgement', '',...
    ...
    'contributor_name', 'Joost den Bieman',...
    'contributor_role', 'project team member',...
    ...
    'publisher_name', 'Kees den Heijer',...
    'publisher_url', 'http://www.deltares.nl',...
    'publisher_email', 'kees.denheijer@deltares.nl',...
    ...
    'processing_level', 'preliminary',...
    'license', 'These data can be used freely for research purposes provided that the following source is acknowledged: DELTARES.',...
    'disclaimer', 'This data is made available in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.',...
    },...
    'Format','classic');

ncwriteschema(ncfile, schema)

%%
[geospatial_vertical_min, geospatial_vertical_max] = deal([]);

ncwrite(ncfile, 'time', tr.time)
ncwrite(ncfile, 'id', id)

multiWaitbar( 'CloseAll' );
for i = 1:n
    multiWaitbar(...
        'volumecalculation', i/n, ...
        'label', sprintf('transect %i', id(i)));
    tr0 = jarkus_transects(...
        'id', id(i),...
        'url', url0,...
        'output', {'cross_shore','min_cross_shore_measurement', 'max_cross_shore_measurement'});
    yearid = ~isnan(tr0.min_cross_shore_measurement);
    bndslimit = nan(sum(yearid),2);
    bndslimit(:,1) = tr0.cross_shore(tr0.min_cross_shore_measurement(yearid));
    bndslimit(:,2) = tr0.cross_shore(tr0.max_cross_shore_measurement(yearid));

    try
        [V, Vc, d, dc, bnds, extrbnds, lon, lat, dlimits] = calc_profilefluct('id', id(i), 'yearid', yearid, 'bndslimit', bndslimit);
        ncwrite(ncfile, 'volume', V, [i 1], [1 1])
        ncwrite(ncfile, 'mean_volume', nanmean(V), i, 1)
        ncwrite(ncfile, 'volumeC', Vc, [i 1], [1 1])
        ncwrite(ncfile, 'mean_volumeC', nanmean(Vc), i, 1)
        ncwrite(ncfile, 'bounds', reshape(bnds, [1,nt,2]), [i 1 1], [1 1 1])
        ncwrite(ncfile, 'boundsC', extrbnds, [i 1], [1 1])
        ncwrite(ncfile, 'boundsC_lon', lon(:)', [i 1], [1 1])
        ncwrite(ncfile, 'boundsC_lat', lat(:)', [i 1], [1 1])
        ncwrite(ncfile, 'contourdepth', d, [i 1], [1 1])
        ncwrite(ncfile, 'contourdepthC', dc, [i 1], [1 1])
        ncwriteatt(ncfile, '/', 'date_modified', datestr(nowutc, datefmt))
        geospatial_vertical_min = min([geospatial_vertical_min; dlimits(:)]);
        geospatial_vertical_max = max([geospatial_vertical_max; dlimits(:)]);
    catch
        fprintf('error in transect %i', id(i));
    end
end

%% complete global attributes
lat_all = ncread(ncfile, 'boundsC_lat');
ncwriteatt(ncfile, '/', 'geospatial_lat_min', min(lat_all(:)));
ncwriteatt(ncfile, '/', 'geospatial_lat_max', max(lat_all(:)));
ncwriteatt(ncfile, '/', 'geospatial_lat_resolution', min(diff(sort(lat_all(:)))));
lon_all = ncread(ncfile, 'boundsC_lon');
ncwriteatt(ncfile, '/', 'geospatial_lon_min', min(lon_all(:)));
ncwriteatt(ncfile, '/', 'geospatial_lon_max', max(lon_all(:)));
ncwriteatt(ncfile, '/', 'geospatial_lon_resolution', min(diff(sort(lon_all(:)))));
ncwriteatt(ncfile, '/', 'geospatial_vertical_min', geospatial_vertical_min);
ncwriteatt(ncfile, '/', 'geospatial_vertical_max', geospatial_vertical_max);
