% Example: create a working 2D plot with time dimension
% 0) Define the variables -> you get this from you model results or data
timerow = "time vector
zstart 	= "wave heights with 3 dimensions: 1:time, 2: lon, 3: lat
x 	= "x-coordinates"
y 	= "y-coordaintes" (both in local coordinate system)

%% 1) Get all data in right coordinate system + orientation
for ii = 1:16

	% Isolate value you would like to present 
	ztmp1 = squeeze(zstart(xx,ii,:,:));

	% Coordinate conversion
	[lon, lat] = convertCoordinates(x,y,'CS1.name','WGS 84 / UTM zone 31N','CS1.type','xy', 'CS2.name', 'WGS 84','CS2.type','geo');
	lonrow = linspace(min(min(lon)), max(max(lon)), 1000);	% lon needs to increase
	latrow = linspace(min(min(lat)), max(max(lat)), 1500);	% lat needs to increase
	[LONROW, LATROW] = meshgrid(lonrow,latrow);

	% Interpolate to new grid
	id = isnan(ztmp1); ztmp1(id) = 0; id = ~isnan(ztmp1); 
	F = scatteredInterpolant(lon(id), lat(id), ztmp1(id), 'linear', 'none');
	ztmp2 = F(LONROW, LATROW);
	z(ii,:,:) = ztmp2;
end

% 2) Get right dimensions for webserver
mylonsize 	= length(lonrow);     mylatsize = length(latrow); 
mytimesize 	= length(timerow)
znew 		= permute(z, [3 2 1]);		! flip z-values, otherwise doesnt work

%% 3) Create netcdf
% Open the file
dir = 'd:\projects\DSD\netcdf\'; cd(dir);
ncid = netcdf.create(['./waveheight_Tillingham_vegetation_no.nc'],'NC_WRITE'); end
globalID = netcdf.getConstant('NC_GLOBAL');

% Add attributes global to the dataset
netcdf.putAtt(ncid,globalID, 'title',           'Wave height as calculated with XBeach-VEG for 1/100 year storm');
netcdf.putAtt(ncid,globalID, 'comment',         'Simulation is a non-hydrostatic XBeach-VEG simulation with and without vegetation. For more information about the possibilities of XBeach-VEG see https://publicwiki.deltares.nl/display/VegMod/XBeach-VEG');
netcdf.putAtt(ncid,globalID, 'institution',     'Deltares');
netcdf.putAtt(ncid,globalID, 'email',           'kees.nederhoff@deltares.nl');
netcdf.putAtt(ncid,globalID, 'scenario',           'current situation without vegetation'); end
netcdf.putAtt(ncid,globalID, 'URL',           'https://www.youtube.com/watch?v=4G4h7Y3rOKY');
netcdf.putAtt(ncid,globalID, 'terms_for_use',   'Use as you like');
netcdf.putAtt(ncid,globalID, 'disclaimer',      'These data are 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');

% Define the dimensions
dimidtime = netcdf.defDim(ncid,'time',mytimesize);
dimidlat = netcdf.defDim(ncid,'lat',mylatsize);
dimidlon = netcdf.defDim(ncid,'lon',mylonsize);
dimidproj = netcdf.defDim(ncid,'projection',1);

% Define IDs for the dimension variables (pressure,time,latitude,...)
time_ID=netcdf.defVar(ncid,'time','double',[dimidtime]);
latitude_ID=netcdf.defVar(ncid,'lat','double',[dimidlat]);
longitude_ID=netcdf.defVar(ncid,'lon','double',[dimidlon]);
project_ID=netcdf.defVar(ncid,'projection','char',[dimidproj]);

%Define the main variable ()
inundationID = netcdf.defVar(ncid,'Hs','double',[dimidlon dimidlat dimidtime]);

%% Standard names
% Standard names - 1
netcdf.putAtt(ncid,time_ID,'units','Days since 1970-01-01 00:00:00');
netcdf.putAtt(ncid,time_ID,'calendar','gregorian');
netcdf.putAtt(ncid,time_ID,'standard_name','time');
netcdf.putAtt(ncid,time_ID,'long_name','time');

% Standard names - 3
netcdf.putAtt(ncid,latitude_ID,'standard_name','latitude');
netcdf.putAtt(ncid,latitude_ID,'long_name','latitude');
netcdf.putAtt(ncid,latitude_ID,'units','degrees_north');
netcdf.putAtt(ncid,latitude_ID,'axis','Y');

% Standard names - 2
netcdf.putAtt(ncid,longitude_ID,'standard_name','longitude');
netcdf.putAtt(ncid,longitude_ID,'long_name','longitude');
netcdf.putAtt(ncid,longitude_ID,'units','degrees_east');
netcdf.putAtt(ncid,longitude_ID,'axis','X');

% Standard names - 4
netcdf.putAtt(ncid,project_ID,'long_name','wgs84');
netcdf.putAtt(ncid,project_ID,'EPSG_code','EPSG:4326');
netcdf.putAtt(ncid,project_ID,'proj4_params','+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs');
netcdf.putAtt(ncid,project_ID,'grid_mapping_name','latitude_longitude');

% Standard names - 5
netcdf.putAtt(ncid,inundationID,'standard_name','Hs');
netcdf.putAtt(ncid,inundationID,'long_name','wave_height');
netcdf.putAtt(ncid,inundationID,'units','meter');
netcdf.putAtt(ncid,inundationID,'_FillValue',-999);
id = znew < 0.04; znew(id) = -999;

%% 4) Finalize the netcdf
%We are done defining the NetCdf
netcdf.endDef(ncid);

%Then store the dimension variables in
netcdf.putVar(ncid,time_ID,timerow);
netcdf.putVar(ncid,longitude_ID,lonrow);
netcdf.putVar(ncid,latitude_ID,latrow);

%Then store my main variable
netcdf.putVar(ncid,inundationID,znew);

%We're done, close the netcdf
netcdf.close(ncid)
fclose('all');