function lidargrid2netcdf(filename, grid)
%LIDARGRID2NETCDF
%
% lidargrid2netcdf(filename, grid)
%
%See also: NETCDF2GRID
%
%   written be Ben deSonniville
%   edited by afarris@usgs.gov, May 2010
%       made small changes to comments and instead of setting size of
%       'time' to be 2, I had set it to the size of grid.time
%       Also made other variables (like shorePos) functions of time.

STRINGSIZE = 100;
%% Create file
mode = bitor ( nc_clobber_mode, nc_64bit_offset_mode );
nc_create_empty ( filename , mode )
% nc_create_empty(filename)
nc_padheader ( filename, 200000 );
%% Put global attributes
nc_attput( filename, nc_global, 'title', 'Lidar Data');
nc_attput( filename, nc_global, 'institution', 'USGS');
nc_attput( filename, nc_global, 'source', 'Laser altimetry');
nc_attput( filename, nc_global, 'history', ['data received from USGS, converted to netcdf on ' date]);
nc_attput( filename, nc_global, 'references', 'This data set was collected with the NASA Airborne Topographic Lidar 3 (ATM3) LIDAR (LIght Detection and Ranging) instrument which was designed and developed by the NASA Observational Science Branch (OSB Code 972) Goddard Space Flight Center, Wallops Island VA. The ATM3 measures the time of flight of a narrow pulse (5 nanoseconds) of green laser light (532 nm) from the aircraft to the ground and back to the aircraft. The laser is pulsed at 5000 Hz and the pulses are reflected from a 22 degree off-nadir scanning mirror rotating at 20 Hz which directs the laser light and the telescope view in an elliptical scan pattern under the aircraft. The telescope collects the reflected laser light which passes through a narrow bandpass filter (to reject background light) and is detected by a photomultiplier tube detector. The ATM3 time of flight measurements are calibrated and combined with GPS data (to provide position) and aircraft inertial gyro data (to provide the aircraft attitude when each laser measurement was recorded) resulting in a calculated latitude, longitude, and elevation for each laser spot in the 250 laser shot elliptical scan under the aircraft. The scan pattern (at 22 degrees off nadir) covers a width of ~ 560 meters from a nominal aircraft altitude of 700 meters. Depending on aircraft ground speed, several consecutive scans may overlap along the aircraft flight track. Normally parallel flight tracks will be flown with ~20% overlap to cover areas wider than a single ATM3 scan width. The resultant data set is normally resampled and binned to create a digital elevation model of the desired area.'); % see http://www.csc.noaa.gov/crs/tcm/ldartdat/metatemplate/atm3_template.html
nc_attput( filename, nc_global, 'Conventions', 'CF-1.4');
% UTM zone 10N; see spatialreference.org using wkt. approach used by gdal.
nc_attput( filename, nc_global, 'spatial_ref', 'PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],AUTHORITY["EPSG","32610"],AXIS["Easting",EAST],AXIS["Northing",NORTH]]');

%% Define dimensions
nc_add_dimension(filename, 'time', length(grid.time));
nc_add_dimension(filename, 'alongshore', length(grid.id));
nc_add_dimension(filename, 'cross_shore', length(grid.crossShoreCoordinate));
nc_add_dimension(filename, 'line', 4);
nc_add_dimension(filename, 'stringsize', STRINGSIZE);

%% Define variables
s.Name      = 'id';
s.Nctype    = nc_int;
s.Dimension = {'alongshore'};
s.Attribute = struct('Name', {'long_name', 'comment'}, 'Value', {'identifier', 'profile number'});
nc_addvar(filename, s);

s.Name      = 'areacode';
s.Nctype    = nc_int;
s.Dimension = {'alongshore'};
s.Attribute = struct('Name', {'long_name', 'comment'}, 'Value', {'area code', 'code for the coastal areas as defined by USGS'});
nc_addvar(filename, s);

s.Name      = 'areaname';
s.Nctype    = nc_char;
s.Dimension = {'alongshore', 'stringsize'};
s.Attribute = struct('Name', {'long_name', 'comment'}, 'Value', {'area name', 'names for the coastal areas as defined by USGS'});
nc_addvar(filename, s);

s.Name      = 'alongshore';
s.Nctype    = nc_double;
s.Dimension = {'alongshore'};
s.Attribute = struct('Name', {'long_name', 'comment'}, 'Value', {'identifier','alongshore point'});
nc_addvar(filename, s);

s.Name      = 'cross_shore';
s.Nctype    = nc_double;
s.Dimension = {'cross_shore'};
s.Attribute = struct('Name',{'long_name', 'comment'}, 'Value', {'identifier', 'cross-shore point'});
nc_addvar(filename, s);

s.Name      = 'contour';
s.Nctype    = nc_double;
s.Dimension = {'alongshore' 'line'}; % Fedor has more universal approach (CF range...)
s.Attribute = struct('Name', {'long_name', 'units', 'comment'}, 'Value', {'outline coordinate', 'm', 'outline coordinate '});
nc_addvar(filename, s);

s.Name      = 'time';
s.Nctype    = nc_double; % nc_int;
s.Dimension = {'time'};
s.Attribute = struct('Name', {'standard_name', 'units', 'comment', 'axis'}, 'Value', {'time', 'year', 'measurement year', 'T'});
nc_addvar(filename, s);


%     [CoordinateSystems, Operations , CoordSysCart ,CoordSysGeo] = GetCoordinateSystems();
%     epsg = 32610; % only input
%     crs = CoordinateSystems([CoordinateSystems.coord_ref_sys_code] == epsg);
%     transform = Operations([Operations.coord_op_code ] == crs.projection_conv_code);
%     s.Name      = 'crs';
%     s.Nctype    = nc_int;
%     s.Dimension = {};
%     s.Attribute = struct('Name', ...
%         {'grid_mapping_name', ...
%         'semi_major_axis', ...
%         'semi_minor_axis', ...
%         'inverse_flattening', ...
%         'latitude_of_projection_origin', ...
%         'longitude_of_projection_origin', ...
%         'false_easting', ...
%         'false_northing', ...
%         'scale_factor_at_projection_origin'}, ...
%         'Value', ...
%         {transform.coordinate_operation_method,...
%         crs.ellipsoid.semi_major_axis, ...
%         crs.ellipsoid.semi_minor_axis, ...
%         crs.ellipsoid.inv_flattening, ...
%         transform.parameters(strcmp({transform.parameters.name},'Latitude_of_natural_origin')).value, ...
%         transform.parameters(strcmp({transform.parameters.name},'Longitude_of_natural_origin')).value, ...
%         transform.parameters(strcmp({transform.parameters.name},'False_easting')).value, ...
%         transform.parameters(strcmp({transform.parameters.name},'False_northing')).value, ...
%         transform.parameters(strcmp({transform.parameters.name},'Scale_factor_at_natural_origin')).value});
%     nc_addvar(filename, s);

s.Name      = 'angle';
s.Nctype    = nc_double;
s.Dimension = {'alongshore'};
s.Attribute = struct('Name', {'long_name', 'units', 'comment'}, 'Value', {'angle of transect', 'mradian', 'positive counter clockwise 0 east'});
nc_addvar(filename, s);

s.Name      = 'mean_high_water';
s.Nctype    = nc_double;
s.Dimension = {'alongshore'};
s.Attribute = struct('Name', {'long_name', 'units', 'comment'}, 'Value', {'mean high water', 'm', 'mean high water'});
nc_addvar(filename, s);

s.Name      = 'mean_low_water';
s.Nctype    = nc_double;
s.Dimension = {'alongshore'};
s.Attribute = struct('Name', {'long_name', 'units', 'comment'}, 'Value', {'mean low water', 'm', 'mean low water'});
nc_addvar(filename, s);

% Lidar specific variables

s.Name      = 'shorepos';
s.Nctype    = nc_double;
s.Dimension = {'alongshore'};
s.Attribute = struct('Name', {'long_name', 'units', 'comment'}, 'Value', {'shoreline position', 'm', 'shoreline position'});
nc_addvar(filename, s);

s.Name      = 'shore_north';
s.Nctype    = nc_double;
s.Dimension = {'alongshore'};
s.Attribute = struct('Name', {'long_name', 'units', 'comment'}, 'Value', {'Northing of shoreline position', 'm', 'Northing of shoreline position'});
nc_addvar(filename, s);

s.Name      = 'shore_east';
s.Nctype    = nc_double;
s.Dimension = {'alongshore'};
s.Attribute = struct('Name', {'long_name', 'units', 'comment'}, 'Value', {'Easting of shoreline position', 'm', 'Easting of shoreline position'});
nc_addvar(filename, s);

s.Name      = 'significant_wave_height';
s.Nctype    = nc_double;
s.Dimension = {'alongshore'};
s.Attribute = struct('Name', {'long_name', 'units', 'comment'}, 'Value', {'Significant waveheight', 'm', 'Significant waveheight'});
nc_addvar(filename, s);

s.Name      = 'deep_water_wave_length';
s.Nctype    = nc_double;
s.Dimension = {'alongshore'};
s.Attribute = struct('Name', {'long_name', 'units', 'comment'}, 'Value', {'Deep water wave length', 'm', 'Deep water wave length'});
nc_addvar(filename, s);

s.Name      = 'beach_slope';
s.Nctype    = nc_double;
s.Dimension = {'alongshore'};
s.Attribute = struct('Name', {'long_name', 'units', 'comment'}, 'Value', {'Beach slope', 'degree', 'Beach slope'});
nc_addvar(filename, s);

s.Name      = 'bias';
s.Nctype    = nc_double;
s.Dimension = {'alongshore'};
s.Attribute = struct('Name', {'long_name', 'units', 'comment'}, 'Value', {'Bias', '-', 'Bias'});
nc_addvar(filename, s);

s.Name      = 'shore_north_1880';
s.Nctype    = nc_double;
s.Dimension = {'alongshore'};
s.Attribute = struct('Name', {'long_name', 'units', 'comment'}, 'Value', {'Northing of shoreline position of 1880', 'm', 'Northing of shoreline position of 1880'});
nc_addvar(filename, s);

s.Name      = 'shore_east_1880';
s.Nctype    = nc_double;
s.Dimension = {'alongshore'};
s.Attribute = struct('Name', {'long_name', 'units', 'comment'}, 'Value', {'Easting of shoreline position of 1880', 'm', 'Easting of shoreline position of 1880'});
nc_addvar(filename, s);

s.Name      = 'shore_north_1930';
s.Nctype    = nc_double;
s.Dimension = {'alongshore'};
s.Attribute = struct('Name', {'long_name', 'units', 'comment'}, 'Value', {'Northing of shoreline position of 1930', 'm', 'Northing of shoreline position of 1930'});
nc_addvar(filename, s);

s.Name      = 'shore_east_1930';
s.Nctype    = nc_double;
s.Dimension = {'alongshore'};
s.Attribute = struct('Name', {'long_name', 'units', 'comment'}, 'Value', {'Easting of shoreline position of 1930', 'm', 'Easting of shoreline position of 1930'});
nc_addvar(filename, s);

%% Store index variables
nc_varput(filename, 'time', grid.time, 0, length(grid.time));
nc_varput(filename, 'id', grid.id');
nc_varput(filename, 'areacode', grid.areaCode');

areanames = grid.areaName;
areanames(:, size(areanames,2)+1:STRINGSIZE) = ' ';

nc_varput(filename, 'areaname', areanames);
nc_varput(filename, 'alongshore', grid.alongshoreCoordinate');
nc_varput(filename, 'cross_shore', grid.crossShoreCoordinate');
nc_varput(filename, 'contour', grid.contour');

% Lidar specific variables

nc_varput(filename, 'mean_high_water', grid.meanHighWater');
nc_varput(filename, 'shorepos', grid.shorepos');
nc_varput(filename, 'shore_east', grid.shoreEast');
nc_varput(filename, 'shore_north', grid.shoreNorth');

nc_varput(filename, 'significant_wave_height', grid.significant_wave_height');
nc_varput(filename, 'deep_water_wave_length', grid.deep_water_wave_length');
nc_varput(filename, 'beach_slope', grid.beach_slope');
nc_varput(filename, 'bias', grid.bias');
%% these are big therefor done seperate

s.Name      = 'altitude';
s.Nctype    = nc_double;
s.Dimension = {'time', 'alongshore', 'cross_shore'};
s.Attribute = struct('Name', {'standard_name', 'units', 'comment', 'coordinates', 'grid_mapping', '_FillValue'}, ...
    'Value', {'surface_altitude', 'm', 'altitude above geoid (NAVD 88)', 'lat lon', 'crs', -9999});
nc_addvar(filename, s);

s.Name      = 'cross_shore_distance';
s.Nctype    = nc_double;
s.Dimension = {'time', 'alongshore', 'cross_shore'};
s.Attribute = struct('Name', {'standard_name', 'units', 'axis'}, 'Value', {'crossshore distance', 'm', 'X'});
nc_addvar(filename, s);

s.Name      = 'lat';
s.Nctype    = nc_double;
s.Dimension = {'time' 'alongshore', 'cross_shore'};
s.Attribute = struct('Name', {'standard_name', 'units', 'axis'}, 'Value', {'latitude', 'degree_north', 'X'});
nc_addvar(filename, s);

s.Name      = 'lon';
s.Nctype    = nc_double;
s.Dimension = {'time' 'alongshore', 'cross_shore'};
s.Attribute = struct('Name', {'standard_name', 'units', 'axis'}, 'Value', {'longitude', 'degree_east', 'Y'});
nc_addvar(filename, s);

s.Name      = 'x';
s.Nctype    = nc_double;
s.Dimension = {'time' 'alongshore', 'cross_shore'};
s.Attribute = struct('Name', {'standard_name', 'units', 'axis'}, 'Value', {'projection_x_coordinate', 'm', 'X'});
nc_addvar(filename, s);

s.Name      = 'y';
s.Nctype    = nc_double;
s.Dimension = {'time' 'alongshore', 'cross_shore'};
s.Attribute = struct('Name', {'standard_name', 'units', 'axis'}, 'Value', {'projection_y_coordinate', 'm', 'Y'});
nc_addvar(filename, s);

s.Name      = 'regression';
s.Nctype    = nc_double;
s.Dimension = {'time' 'alongshore', 'cross_shore'};
s.Attribute = struct('Name', {'standard_name', 'units', 'axis'}, 'Value', {'Regression flag', 'm', 'Regression value'});
nc_addvar(filename, s);


end % function grid2netcdf

