%% Create netCDF-CF of curvilinear lat-lon grid
%
% example of how to make a netCDF file with CF conventions of a
% variable that is defined on a grid that is curvi-linear
% in a lat-lon coordinate system. In this case
% the dimensions (m,n) do not coincide with the coordinate axes.
%
% This case is described in:
% http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html#id2984605
% as "Two-Dimensional Latitude, Longitude, Coordinate Variables".
%
% An example of a curvi-linear lat,lon grid is for instance a satellite
% image: a digital image with a constant width in terms of number of
% pixels ncols (swath), and a length nrows determined by how long a ground
% station can receive the satellite above the (electro-magnetic) horizon.
%
% ^ latitude (degrees_north)
% |
% | / /\
% |ncols (Track)/ /15\ \
% | / / \ \
% | / /10 \ \
% | / / \ \
% | (5 14\ \
% | \ 9 \ \
% | \4 \ |
% | \ \ |
% | )3 8 13) | nrows (Xtrack)
% | / / |
% | /2 / |
% | / 7 / /
% | (1 12/ /
% | \ / /
% | \6 / /
% | \ / /
% | \11/ /
% | \/
% +----------------------> longitude
% (degrees_east)
%
% In MATLAB, a 2D variable varies fastest along the rows, so we will borrow
% from HDF-EOS swath terminology and call that dimension 'Xtrack'. The
% along track dimension would then be along the columns and we will call
% this 'Track'. This will determine our order of dimensions. Note that if
% we use row-major order conventions, the order of dimensions will be
% reversed.
%
% Note that ncBrowse does not contain plot support for
% curvi-linear grids, so ncBrowse will display the same
% rectangular plot as for the netCDF file created by
% NC_CF_GRID_WRITE_LAT_LON_ORTHOGONAL_TUTORIAL, albeit with
% different axes annotations (col/row instead of lat/lon).
%
%See also: SNCTOOLS, NC_CF_GRID_WRITE_LAT_LON_ORTHOGONAL_TUTORIAL,
%% Define meta-info: global
OPT.title = 'Example of curvilinear grid conforming to CF conventions.';
OPT.institution = 'TU Delpht';
OPT.source = 'example data';
OPT.history = ['tranformation to netCDF: $HeadURL: https://repos.deltares.nl/repos/OpenEarthTools/trunk/matlab/io/netcdf/nctools/nc_cf_grid_write_lat_lon_curvilinear_tutorial.m $'];
OPT.references = 'http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html';
OPT.email = 'john.g.evans.ne@gmail.com';
OPT.comment = 'Adapted from example provided by Gerben de Boer.';
OPT.version = '1.0';
OPT.acknowledge =['These data can be used freely for research purposes provided that the following source is acknowledged: ',OPT.institution];
OPT.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.';
%% Define dimensions/coordinates: lat,lon matrices
lon1 = [2 4 6];
lat1 = [50 51 52 53 54];
[lat2,lon2] = ndgrid(lat1,lon1);
ang = [-15 -15 -15; 0 0 0; 15 15 15; 30 30 30; 45 45 45];
OPT.lat = lat2 + sind(ang).*lon2./2;
OPT.lon = lon2 + cosd(ang).*lon2./2; clear lon1 lon2 lat1 lat2
OPT.ncols = size(OPT.lon,2);
OPT.nrows = size(OPT.lat,1);
OPT.lat_type = 'single'; % 'single', 'double' for high-resolution data (eps 1m)
OPT.lon_type = 'single'; % 'single', 'double' for high-resolution data (eps 1m)
OPT.wgs84.code = 4326; % epsg code of global grid: http://www.epsg-registry.org/
OPT.wgs84.name = 'WGS 84';
OPT.wgs84.semi_major_axis = 6378137.0;
OPT.wgs84.semi_minor_axis = 6356752.314247833;
OPT.wgs84.inv_flattening = 298.2572236;
%% Define variable (define some data)
OPT.val = [1 2 3 4 5;6 7 8 9 10;11 12 nan 14 15]'; % use ncols as 1st array dimension to get correct plot in ncBrowse (snctools swaps for us)
OPT.varname = 'depth'; % free to choose: will appear in netCDF tree
OPT.units = 'm'; % from UDunits package: http://www.unidata.ucar.edu/software/udunits/
OPT.long_name = 'bottom depth';% free to choose: will appear in plots
OPT.standard_name = 'sea_floor_depth_below_geoid'; % or 'altitude'
OPT.val_type = 'single'; % 'single' or 'double'
OPT.fillvalue = single(-9999);
%% 1.a Create netCDF file
ncfile = fullfile(fileparts(mfilename('fullpath')),[mfilename,'.nc']);
nc_create_empty (ncfile)
%% 1.b Add overall meta info
% http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html#description-of-file-contents
nc_attput(ncfile, nc_global, 'title' , OPT.title);
nc_attput(ncfile, nc_global, 'institution' , OPT.institution);
nc_attput(ncfile, nc_global, 'source' , OPT.source);
nc_attput(ncfile, nc_global, 'history' , OPT.history);
nc_attput(ncfile, nc_global, 'references' , OPT.references);
nc_attput(ncfile, nc_global, 'email' , OPT.email);
nc_attput(ncfile, nc_global, 'comment' , OPT.comment);
nc_attput(ncfile, nc_global, 'version' , OPT.version);
nc_attput(ncfile, nc_global, 'Conventions' , 'CF-1.4');
nc_attput(ncfile, nc_global, 'CF:featureType', 'Grid'); % https://cf-pcmdi.llnl.gov/trac/wiki/PointObservationConventions
nc_attput(ncfile, nc_global, 'terms_for_use' , OPT.acknowledge);
nc_attput(ncfile, nc_global, 'disclaimer' , OPT.disclaimer);
%% 2 Create matrix span dimensions
% http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html#dimensions
nc_adddim(ncfile, 'Xtrack', OPT.nrows); % !!! use this as 1st array dimension to get correct plot in ncBrowse (snctools swaps for us)
nc_adddim(ncfile, 'Track', OPT.ncols); % !!! use this as 2nd array dimension to get correct plot in ncBrowse (snctools swaps for us)
% You might insert a vector 'col' that runs [OPT.ncols:-1:1] to have
% the arcGIS ASCII file approach of having upper-left corner of
% the data matrix at index (1,1) rather than the default of having the
% lower-left corner of the data matrix at index (1,1).
% nc_add_dimension(ncfile, 'time', 1); % if you would like to include more instances of the same grid,
% you can optionally use 'time' as a 3rd dimension.
% T,(Z),Y,X is the recommended dimension CF order.
%% 3.a Create coordinate variables: longitude
% http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html#longitude-coordinate
clear nc;ifld = 1;
nc(ifld).Name = 'lon';
nc(ifld).Datatype = OPT.lon_type;
nc(ifld).Dimension = {'Xtrack','Track'};
nc(ifld).Attribute( 1) = struct('Name', 'long_name' ,'Value', 'longitude');
nc(ifld).Attribute(end+1) = struct('Name', 'units' ,'Value', 'degrees_east');
nc(ifld).Attribute(end+1) = struct('Name', 'standard_name' ,'Value', 'longitude');
nc(ifld).Attribute(end+1) = struct('Name', 'actual_range' ,'Value', [min(OPT.lon(:)) max(OPT.lon(:))]);
nc(ifld).Attribute(end+1) = struct('Name', 'coordinates' ,'Value', 'lat lon'); % !!! lon matrix can be plotted as a function of lat and itself
nc(ifld).Attribute(end+1) = struct('Name', 'grid_mapping' ,'Value', 'wgs84');
%% 3.b Create coordinate variables: latitude
% http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html#latitude-coordinate
ifld = ifld + 1;
nc(ifld).Name = 'lat';
nc(ifld).Datatype = OPT.lat_type;
nc(ifld).Dimension = {'Xtrack','Track'};
nc(ifld).Attribute( 1) = struct('Name', 'long_name' ,'Value', 'latitude');
nc(ifld).Attribute(end+1) = struct('Name', 'units' ,'Value', 'degrees_north');
nc(ifld).Attribute(end+1) = struct('Name', 'standard_name' ,'Value', 'latitude');
nc(ifld).Attribute(end+1) = struct('Name', 'actual_range' ,'Value', [min(OPT.lat(:)) max(OPT.lat(:))]);
nc(ifld).Attribute(end+1) = struct('Name', 'coordinates' ,'Value', 'lat lon'); % !!! lat matrix can be plotted as a function of lon and itself
nc(ifld).Attribute(end+1) = struct('Name', 'grid_mapping' ,'Value', 'wgs84');
%% 3.c Create coordinate variables: coordinate system: WGS84 default
% global ellispes: WGS 84, ED 50, INT 1924, ETRS 89 and the upcoming ETRS update etc.
% http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html#grid-mappings-and-projections
% http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html#appendix-grid-mappings
ifld = ifld + 1;
nc(ifld).Name = 'wgs84'; % preferred
nc(ifld).Datatype = nc_int;
nc(ifld).Dimension = {};
nc(ifld).Attribute = struct('Name', ...
{'name',...
'epsg',...
'grid_mapping_name',...
'semi_major_axis', ...
'semi_minor_axis', ...
'inverse_flattening', ...
'comment'}, ...
'Value', ...
{OPT.wgs84.name,...
OPT.wgs84.code,...
'latitude_longitude',...
OPT.wgs84.semi_major_axis, ...
OPT.wgs84.semi_minor_axis, ...
OPT.wgs84.inv_flattening, ...
'value is equal to EPSG code'});
%% 4 Create dependent variable
% http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html#variables
% Parameters with standard names:
% http://cf-pcmdi.llnl.gov/documents/cf-standard-names/standard-name-table/current/
ifld = ifld + 1;
nc(ifld).Name = OPT.varname;
nc(ifld).Datatype = OPT.val_type;
nc(ifld).Dimension = {'Xtrack','Track'}; % {'time','col','row'}
nc(ifld).Attribute( 1) = struct('Name', 'long_name' ,'Value', OPT.long_name );
nc(ifld).Attribute(end+1) = struct('Name', 'units' ,'Value', OPT.units );
nc(ifld).Attribute(end+1) = struct('Name', '_FillValue' ,'Value', OPT.fillvalue );
nc(ifld).Attribute(end+1) = struct('Name', 'actual_range' ,'Value', [min(OPT.val(:)) max(OPT.val(:))]);
nc(ifld).Attribute(end+1) = struct('Name', 'coordinates' ,'Value', 'lat lon');
nc(ifld).Attribute(end+1) = struct('Name', 'grid_mapping' ,'Value', 'wgs84');
if ~isempty(OPT.standard_name)
nc(ifld).Attribute(end+1) = struct('Name', 'standard_name' ,'Value', OPT.standard_name);
end
%% 5.a Create all variables with attributes
for ifld=1:length(nc)
nc_addvar(ncfile, nc(ifld));
end
%% 5.b Fill all variables
nc_varput(ncfile, 'lon' , OPT.lon );
nc_varput(ncfile, 'lat' , OPT.lat );
nc_varput(ncfile, 'wgs84' , OPT.wgs84.code);
nc_varput(ncfile, OPT.varname , OPT.val );
%% 6 Check file summary
nc_dump(ncfile);
fid = fopen(fullfile(fileparts(mfilename('fullpath')),[mfilename,'.cdl']),'w');
nc_dump(ncfile,fid);
fclose(fid)
%% 7 Load the data: using the variable names from nc_dump
Da.dep = nc_varget(ncfile,'depth');
Da.lat = nc_varget(ncfile,'lon');
Da.lon = nc_varget(ncfile,'lat')
Running this example will produce output that looks as follows:
netCDF /Users/jevans/Downloads/nc_cf_grid_write_lat_lon_curvilinear_tutorial.nc {
dimensions:
Xtrack = 5 ;
Track = 3 ;
variables:
// Preference 'PRESERVE_FVD': true,
// dimensions consistent with native MATLAB netcdf package, not with ncBrowse.
single lon(Xtrack,Track), shape = [5 3]
lon:long_name = "longitude"
lon:units = "degrees_east"
lon:standard_name = "longitude"
lon:actual_range = 2.70711 9
lon:coordinates = "lat lon"
lon:grid_mapping = "wgs84"
single lat(Xtrack,Track), shape = [5 3]
lat:long_name = "latitude"
lat:units = "degrees_north"
lat:standard_name = "latitude"
lat:actual_range = 49.2235 56.1213
lat:coordinates = "lat lon"
lat:grid_mapping = "wgs84"
int32 wgs84([]), shape = [1]
wgs84:name = "WGS 84"
wgs84:epsg = 4326
wgs84:grid_mapping_name = "latitude_longitude"
wgs84:semi_major_axis = 6.37814e+06
wgs84:semi_minor_axis = 6.35675e+06
wgs84:inverse_flattening = 298.257
wgs84:comment = "value is equal to EPSG code"
single depth(Xtrack,Track), shape = [5 3]
depth:long_name = "bottom depth"
depth:units = "m"
depth:_FillValue = -9999.000000 f
depth:actual_range = 1 15
depth:coordinates = "lat lon"
depth:grid_mapping = "epsg"
depth:standard_name = "sea_floor_depth_below_geoid"
//global attributes:
:title = "Example of curvilinear grid conforming to CF conventions."
:institution = "TU Delpht"
:source = "example data"
:history = "tranformation to netCDF: $HeadURL: https://repos.deltares.nl/repos/OpenEarthTools/trunk/matlab/io/netcdf/nctools/nc_cf_grid_write_lat_lon_curvilinear_tutorial.m $"
:references = "http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html"
:email = "john.g.evans.ne@gmail.com"
:comment = "Adapted from example provided by Gerben de Boer."
:version = "1.0"
:Conventions = "CF-1.4"
:CF:featureType = "Grid"
:terms_for_use = "These data can be used freely for research purposes provided that the following source is acknowledged: TU Delpht"
: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."
}