Contents

Create netCDF-CF of orthogonal x-y grid

example of how to make a netCDF file with CF conventions of a
variable that is defined on a grid that is orthogonal
in a x-y coordinate system. In this special case
the dimensions coincide with the x-y coordinate axes.
The associated lat-lon coordinate matrices are curvi-linear.
This case is described in:
http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html#grid-mappings-and-projections
as "Horizontal Coordinate Reference Systems, Grid Mappings, and Projections".
An example of an orthogonal x,y grid is for instance
a data product defined on an otrhogonal grid in a certain
local geograhic projection by programs like arcGIS. For each
vertex in this grid a the lat,lon values need to be calculated
to project it on the globe (e.g. Google Earth).
  ^ latitude (degrees_north)           ^ y(m)
|         x                          |
| ncols /   \                        |        ncols
|     /  /\   \              coordinate	  +------+
|   /   /15\    \          transformation
|  x  /10   \     \       <==============>	 +--------+
|    <5     14\     \                |      |05 10 15|  +
|     \   9    \      \              |      |        |  |
|      \4       \      |             |      |04 09 14|  |
|       \        \     |             |      |        |  |
|        )3  8  xx)    | nrows       |      |03 08 xx|  | nrows
|       /        /     |             |      |        |  |
|      /2       /      |             |      |02 07 12|  |
|     /   7    /      /              |      |        |  |
|    <1     12/     /                |      |01 06 11|  +
|     \6    /     /                  |      +--------+
|       \11/    /                    |
|        \/   x                      |
|                                    |
+----------------------> longitude   +----------------------> x
(degrees_east)                          (m)

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, NC_CF_GRID_WRITE,
%          NC_CF_GRID_WRITE_LAT_LON_ORTHOGONAL_TUTORIAL,
%          NC_CF_GRID_WRITE_LAT_LON_CURVILINEAR_TUTORIAL,
%          NC_CF_GRID_WRITE_X_Y_CURVILINEAR_TUTORIAL,
%          NC_CF_STATIONTIMESERIES

% This tool is part of <a href="http://www.OpenEarth.eu">OpenEarthTools</a> under the <a href="http://www.gnu.org/licenses/gpl.html">GPL</a> license.

Define meta-info: global: x,y sticks > lat,lon matrices

   OPT.title                  = '';
OPT.institution            = '';
OPT.source                 = '';
OPT.history                = ['tranformation to netCDF: $HeadURL: https://repos.deltares.nl/repos/OpenEarthTools/trunk/matlab/io/netcdf/nctools/nc_cf_grid_write_x_y_orthogonal_tutorial.m $'];
OPT.references             = '';
OPT.email                  = '';
OPT.comment                = '';
OPT.version                = '';
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

   OPT.x                      = [400 500 600].*1e3;
OPT.y                      = [5500:100:5900].*1e3;

[x,y]                       = ndgrid(OPT.x,OPT.y);

OPT.ncols                  = length(OPT.x);
OPT.nrows                  = length(OPT.y);
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.epsg.code              = 32631; % epsg code of local projection
OPT.wgs84.code             = 4326;  % epsg code of global grid
% http://www.epsg-registry.org/
% in the case of a grid defined in a local x-y
% projection, the properties of the grid in a WGS84
% lat,lon system do not have to be specified here, but
% can be retrieved from the log of the coordinate
% transformation carried out by convertCoordinates.
% get (lat,lon) associated with each vertex (x,y), note order [OPT.lon,OPT.lat ...

[OPT.lon,OPT.lat,log]       = convertCoordinates(x,y,'CS1.code',OPT.epsg.code,'CS2.code',OPT.wgs84.code);

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              = nan;

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.a Create matrix span dimensions

    http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html#dimensions
   nc_add_dimension(ncfile, 'x', OPT.ncols); % use this as 1st array dimension to get correct plot in ncBrowse (snctools swaps for us)
nc_add_dimension(ncfile, 'y', OPT.nrows); % use this as 2nd array dimension to get correct plot in ncBrowse (snctools swaps for us)

% You might insert a vector 'col' that runs max(x):-dx:min(x) 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. see
% nc_cf_stationTimeSeries_write_tutorial for info on time.

3.a Create coordinate vectors: x and y

   clear nc;ifld = 1;
nc(ifld).Name             = 'x'; % dimension 'x' is here filled with variable 'x'
nc(ifld).Nctype           = nc_type(OPT.lon_type);
nc(ifld).Dimension        = {'x'};
nc(ifld).Attribute(    1) = struct('Name', 'long_name'      ,'Value', 'x Rijksdriehoek');
nc(ifld).Attribute(end+1) = struct('Name', 'units'          ,'Value', 'm');
nc(ifld).Attribute(end+1) = struct('Name', 'standard_name'  ,'Value', 'projection_x_coordinate'); % standard name
nc(ifld).Attribute(end+1) = struct('Name', 'actual_range'   ,'Value', [min(OPT.x(:)) max(OPT.x(:))]);
nc(ifld).Attribute(end+1) = struct('Name', 'grid_mapping'   ,'Value', 'epsg');

ifld = ifld + 1;
nc(ifld).Name             = 'y'; % dimension 'y' is here filled with variable 'y'
nc(ifld).Nctype           = nc_type(OPT.lat_type);
nc(ifld).Dimension        = {'y'};
nc(ifld).Attribute(    1) = struct('Name', 'long_name'      ,'Value', 'y Rijksdriehoek');
nc(ifld).Attribute(end+1) = struct('Name', 'units'          ,'Value', 'm');
nc(ifld).Attribute(end+1) = struct('Name', 'standard_name'  ,'Value', 'projection_y_coordinate'); % standard name
nc(ifld).Attribute(end+1) = struct('Name', 'actual_range'   ,'Value', [min(OPT.y(:)) max(OPT.y(:))]);
nc(ifld).Attribute(end+1) = struct('Name', 'grid_mapping'   ,'Value', 'epsg');

3.b Create coordinate variables: coordinate system: epsg

    http://www.epsg-registry.org/
   ifld = ifld + 1;
nc(ifld).Name         = 'epsg';
nc(ifld).Nctype       = nc_int;
nc(ifld).Dimension    = {};
nc(ifld).Attribute    = nc_cf_grid_mapping(OPT.epsg.code);
var2evalstr(nc(ifld).Attribute)

OPT = 

debug: 0
grid_mapping_name: 'transverse_mercator'


ans =

variable1(1).Name = 'name';
variable1(1).Value = 'WGS 84 / UTM zone 31N';
variable1(2).Name = 'epsg';
variable1(2).Value = 32631;
variable1(3).Name = 'epsg_name';
variable1(3).Value = 'Transverse Mercator';
variable1(4).Name = 'grid_mapping_name';
variable1(4).Value = 'transverse_mercator';
variable1(5).Name = 'semi_major_axis';
variable1(5).Value = 6378137;
variable1(6).Name = 'semi_minor_axis';
variable1(6).Value = 6356752.3142;
variable1(7).Name = 'inverse_flattening';
variable1(7).Value = 298.2572;
variable1(8).Name = 'latitude_of_projection_origin';
variable1(8).Value = 0;
variable1(9).Name = 'longitude_of_projection_origin';
variable1(9).Value = 3;
variable1(10).Name = 'false_easting';
variable1(10).Value = 500000;
variable1(11).Name = 'false_northing';
variable1(11).Value = 0;
variable1(12).Name = 'scale_factor_at_projection_origin';
variable1(12).Value = 0.9996;
variable1(13).Name = 'comment';
variable1(13).Value = 'value is equal to EPSG code';

3.c Create coordinate variables: longitude

    http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html#longitude-coordinate
   ifld = ifld + 1;
nc(ifld).Name             = 'lon';
nc(ifld).Nctype           = nc_type(OPT.lon_type);
nc(ifld).Dimension        = {'x','y'};
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');
nc(ifld).Attribute(end+1) = struct('Name', 'grid_mapping'   ,'Value', 'wgs84');

3.d 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).Nctype           = nc_type(OPT.lat_type);
nc(ifld).Dimension        = {'x','y'};
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');
nc(ifld).Attribute(end+1) = struct('Name', 'grid_mapping'   ,'Value', 'wgs84');

3.e 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).Nctype       = nc_int;
nc(ifld).Dimension    = {};
nc(ifld).Attribute    = nc_cf_grid_mapping(OPT.wgs84.code);
var2evalstr(nc(ifld).Attribute)

ans =

variable1(1).Name = 'name';
variable1(1).Value = 'WGS 84';
variable1(2).Name = 'epsg';
variable1(2).Value = 4326;
variable1(3).Name = 'grid_mapping_name';
variable1(3).Value = 'latitude_longitude';
variable1(4).Name = 'semi_major_axis';
variable1(4).Value = 6378137;
variable1(5).Name = 'semi_minor_axis';
variable1(5).Value = 6356752.3142;
variable1(6).Name = 'inverse_flattening';
variable1(6).Value = 298.2572;
variable1(7).Name = 'comment';
variable1(7).Value = '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).Nctype           = nc_type(OPT.val_type);
nc(ifld).Dimension        = {'x','y'}; % {'time','x','y'};
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', 'epsg');
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, 'x'            , OPT.x         );
nc_varput(ncfile, 'y'            , OPT.y         );
nc_varput(ncfile, 'epsg'         , OPT.epsg.code );
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)

netCDF F:\checkouts\OpenEarthTools\matlab\io\netcdf\nctools\nc_cf_grid_write_x_y_orthogonal_tutorial.nc { 

dimensions:
x = 3 ;
y = 5 ;


variables:
// Preference 'PRESERVE_FVD':  false,
// dimensions consistent with ncBrowse, not with native MATLAB netcdf package.
single x(x), shape = [3]
x:long_name = "x Rijksdriehoek" 
x:units = "m" 
x:standard_name = "projection_x_coordinate" 
x:actual_range = 400000 600000 
x:grid_mapping = "epsg" 
single y(y), shape = [5]
y:long_name = "y Rijksdriehoek" 
y:units = "m" 
y:standard_name = "projection_y_coordinate" 
y:actual_range = 5.5e+006 5.9e+006 
y:grid_mapping = "epsg" 
int32 epsg([]), shape = [1]
epsg:name = "WGS 84 / UTM zone 31N" 
epsg:epsg = 32631 
epsg:epsg_name = "Transverse Mercator" 
epsg:grid_mapping_name = "transverse_mercator" 
epsg:semi_major_axis = 6.37814e+006 
epsg:semi_minor_axis = 6.35675e+006 
epsg:inverse_flattening = 298.257 
epsg:latitude_of_projection_origin = 0 
epsg:longitude_of_projection_origin = 3 
epsg:false_easting = 500000 
epsg:false_northing = 0 
epsg:scale_factor_at_projection_origin = 0.9996 
epsg:comment = "value is equal to EPSG code" 
single lon(x,y), shape = [3 5]
lon:long_name = "longitude" 
lon:units = "degrees_east" 
lon:standard_name = "longitude" 
lon:actual_range = 1.50155 4.49845 
lon:coordinates = "lat lon" 
lon:grid_mapping = "wgs84" 
single lat(x,y), shape = [3 5]
lat:long_name = "latitude" 
lat:units = "degrees_north" 
lat:standard_name = "latitude" 
lat:actual_range = 49.6443 53.2493 
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+006 
wgs84:semi_minor_axis = 6.35675e+006 
wgs84:inverse_flattening = 298.257 
wgs84:comment = "value is equal to EPSG code" 
single depth(x,y), shape = [3 5]
depth:long_name = "bottom depth" 
depth:units = "m" 
depth:_FillValue = NaN 
depth:actual_range = 1 15 
depth:coordinates = "lat lon" 
depth:grid_mapping = "epsg" 
depth:standard_name = "sea_floor_depth_below_geoid" 


//global attributes:
:title = "" 
:institution = "" 
:source = "" 
:history = "tranformation to netCDF: $HeadURL: https://repos.deltares.nl/repos/OpenEarthTools/trunk/matlab/io/netcdf/nctools/nc_cf_grid_write_x_y_orthogonal_tutorial.m $" 
:references = "" 
:email = "" 
:comment = "" 
:version = "" 
: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: " 
: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." 
}

ans =

0

7.a 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')

Da = 

dep: [3x5 double]
lat: [3x5 single]
lon: [3x5 single]
datenum: [6x1 double]
var: [1 2 3 NaN 5 6]

7.b Load the data: using standard_names and coordinate attribute

   depname  = nc_varfind(ncfile,'attributename', 'standard_name', 'attributevalue', 'sea_floor_depth_below_geoid')
Db.z     = nc_varget(ncfile,depname);

coords   = nc_attget(ncfile,depname,'coordinates');
[ax1,coords] = strtok(coords); ax2 = strtok(coords);
if strcmpi(nc_attget(ncfile,ax1,'standard_name'),'latitude')
Db.lat   = nc_varget(ncfile,ax1);
Db.lon   = nc_varget(ncfile,ax2)
else
Db.lat   = nc_varget(ncfile,ax2);
Db.lon   = nc_varget(ncfile,ax1)
end

depname =

depth


Db = 

z: [3x5 double]
lat: [3x5 single]
lon: [3x5 single]
datenum: [6x1 double]
var: [1 2 3 NaN 5 6]

7.c Load the data: using a dedicated function developed for grids

   [Dc,Mc] = nc_cf_grid(ncfile,OPT.varname)

Dc = 

lon: [3x5 single]
lat: [3x5 single]
depth: [3x5 double]


Mc = 

lon: [1x1 struct]
lat: [1x1 struct]
depth: [1x1 struct]


this tutorial is based on: nc_cf_grid_write_x_y_orthogonal_tutorial.m (revision: NaN)
Published with MATLAB® 7.7