%KNMI_NOAAPC2NC rewrite binary KNMI NOAA POES SST data files (NOAAPC format) into NetCDF files
%
% Grid data definition, see example:
% http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html#id2984605
%
% The KNMI set contains data from the POES satellites 10 to 18:
% www.knmi.nl
% www.knmi.nl/onderzk/applied/sd/en/AVHRR_archive_KNMI.html
%
%See also: KNMI_NOAAPC_READ
% This tool is part of OpenEarthTools.
% OpenEarthTools is an online collaboration to share and manage data and
% programming tools in an open source, version controlled environment.
% Sign up to recieve regular updates of this function, and to contribute
% your own tools.
%% Version
% $Id$
% $Date$
% $Author$
% $Revision$
% $HeadURL$
% $Keywords: $
% TO DO: find out polar_stereographic parameters to EPSG
% TO DO: save data with scale_factor and add_offset to UINT8 to sace space
% TO DO: add composite info (to cell_methods)
% function knmi_noaapc2nc(...)
%% Initialize
OPT.fillvalue = nan; % NaNs do work in netcdf API
OPT.dump = 0;
OPT.pause = 0;
OPT.debug = 0;
OPT.pack = 0;
OPT.ll = 1;
OPT.refdatenum = datenum(0000,0,0); % matlab datenumber convention: A serial date number of 1 corresponds to Jan-1-0000. Gives wring date sin ncbrowse due to different calenders. Must use doubles here.
OPT.refdatenum = datenum(1970,1,1); % lunix datenumber convention
OPT.timezone = timezone_code2iso('GMT');
%% File loop
OPT.directory.raw = ['F:\checkouts\OpenEarthRawData\knmi\NOAA\noaapc\1990_mom\5\'];
OPT.directory.nc = ['F:\checkouts\OpenEarthRawData\knmi\NOAA\noaapc\mom.nc\1990_mom\5\'];
mkpath(OPT.directory.nc)
OPT.files = dir([OPT.directory.raw filesep '*.SST']);
for ifile=1:length(OPT.files)
OPT.filename = ([OPT.directory.raw, filesep, OPT.files(ifile).name]); % id1-AMRGBVN-196101010000-200801010000.txt
disp(['Processing ',num2str(ifile),'/',num2str(length(OPT.files)),': ',filename(OPT.filename)])
%% 0 Read raw data
D = knmi_noaapc_read(OPT.filename,'center',1,'landmask',nan,'cloudmask',-Inf,'count',OPT.pack); % make sure to set valid_min to prevent -Inf from corrupting color scale in ncBrowse.
D.version = '';
if OPT.debug
pcolorcorcen(D.loncor,D.latcor,D.data)
end
%% 1a Create file
if OPT.pack
OPT.ext = ['_sst_pack',num2str(OPT.pack)];
else
OPT.ext = '';
end
%outputfile = [OPT.directory.nc,filesep, filename(OPT.filename) ,OPT.ext,'.nc']; % 8.3 DOS KNMI name
outputfile = [OPT.directory.nc,filesep,'N',datestr(D.datenum,30),'_SST',OPT.ext,'.nc']; % 30 (ISO 8601) 'yyyymmddTHHMMSS' name
nc_create_empty (outputfile)
%% Add overall meta info
% http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html#description-of-file-contents
%------------------
nc_attput(outputfile, nc_global, 'title' , '');
nc_attput(outputfile, nc_global, 'institution' , 'KNMI');
nc_attput(outputfile, nc_global, 'source' , 'surface observation');
nc_attput(outputfile, nc_global, 'history' ,['Original filename: ',filename(OPT.filename),...
', version:' ,D.version,...
', filedate:',OPT.files(ifile).date,...
', tranformation to NetCDF: $HeadURL$ $Date$ $Author$']);
nc_attput(outputfile, nc_global, 'references' , ',,,');
nc_attput(outputfile, nc_global, 'email' , '');
nc_attput(outputfile, nc_global, 'comment' , '');
nc_attput(outputfile, nc_global, 'version' , D.version);
nc_attput(outputfile, nc_global, 'Conventions' , 'CF-1.4');
nc_attput(outputfile, nc_global, 'CF:featureType' , 'Grid'); % https://cf-pcmdi.llnl.gov/trac/wiki/PointObservationConventions
nc_attput(outputfile, nc_global, 'terms_for_use' , 'These data can be used freely for research purposes provided that the following source is acknowledged: KNMI.');
nc_attput(outputfile, nc_global, '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.');
nc_attput(outputfile, nc_global, 'satnum' , D.satnum);
nc_attput(outputfile, nc_global, 'orbnum' , D.orbnum);
nc_attput(outputfile, nc_global, 'type' , D.type);
nc_attput(outputfile, nc_global, 'yearday' , D.yearday);
%% Add discovery information (test):
% http://www.unidata.ucar.edu/projects/THREDDS/tech/catalog/InvCatalogSpec.html
nc_attput(outputfile, nc_global, 'geospatial_lat_min' , min(D.latcor(:)));
nc_attput(outputfile, nc_global, 'geospatial_lat_max' , max(D.latcor(:)));
nc_attput(outputfile, nc_global, 'geospatial_lon_min' , min(D.loncor(:)));
nc_attput(outputfile, nc_global, 'geospatial_lon_max' , max(D.loncor(:)));
nc_attput(outputfile, nc_global, 'time_coverage_start' , datestr(D.datenum( 1),'yyyy-mm-ddPHH:MM:SS'));
nc_attput(outputfile, nc_global, 'time_coverage_end' , datestr(D.datenum(end),'yyyy-mm-ddPHH:MM:SS'));
nc_attput(outputfile, nc_global, 'geospatial_lat_units' , 'degrees_north');
nc_attput(outputfile, nc_global, 'geospatial_lon_units' , 'degrees_east' );
%% 2 Create dimensions
nc_add_dimension(outputfile, 'time' , 1)
nc_add_dimension(outputfile, 'x' , D.nx)
nc_add_dimension(outputfile, 'y' , D.ny)
nc_add_dimension(outputfile, 'x_cor', D.nx+1)
nc_add_dimension(outputfile, 'y_cor', D.ny+1)
%% 3 Create variables
clear nc
ifld = 0;
%% Coordinate system
% http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html#appendix-grid-mappings
ifld = ifld + 1;
nc(ifld).Name = 'polar_stereographic';
nc(ifld).Nctype = 'char';
nc(ifld).Dimension = {}; % no dimension, dummy variable
nc(ifld).Attribute(end+1) = struct('Name', 'projection_name' ,'Value', 'Polar Stereographic';
nc(ifld).Attribute(end+1) = struct('Name', 'EPSG_code' ,'Value', 'UNKNOWN'); % or UNDEFINED
nc(ifld).Attribute(end+1) = struct('Name', 'proj4_params' ,'Value', D.proj4_params);
% See ADADUC manual for polar_stereographic example
nc(ifld).Attribute( 1) = struct('Name', 'grid_mapping_name' ,'Value', 'polar_stereographic');
nc(ifld).Attribute(end+1) = struct('Name', 'latitude_of_projection_origin' ,'Value', '+90'); % Either +90. or -90.
nc(ifld).Attribute(end+1) = struct('Name', 'straight_vertical_longitude_from_pole','Value', 0);
nc(ifld).Attribute(end+1) = struct('Name', 'scale_factor_at_projection_origin' ,'Value', D.scale_in_m);
nc(ifld).Attribute(end+1) = struct('Name', 'false_easting' ,'Value', 0);
nc(ifld).Attribute(end+1) = struct('Name', 'false_northing' ,'Value', 0);
nc(ifld).Attribute(end+1) = struct('Name', 'inverse_flattening' ,'Value', 298.183263207106);
nc(ifld).Attribute(end+1) = struct('Name', 'semi_major_axis' ,'Value', 6378140000);
nc(ifld).Attribute(end+1) = struct('Name', 'semi_minor_axis' ,'Value', 6356750000);
nc(ifld).Attribute(end+1) = struct('Name', 'longitude_of_prima_meridian' ,'Value', 0);
%% Local Cartesian coordinates
ifld = ifld + 1;
nc(ifld).Name = 'x';
nc(ifld).Nctype = 'int';
nc(ifld).Dimension = {'x'};
nc(ifld).Attribute(1) = struct('Name', 'long_name' ,'Value', 'x-coordinate in Cartesian system');
nc(ifld).Attribute(2) = struct('Name', 'units' ,'Value', 'km');
nc(ifld).Attribute(3) = struct('Name', 'standard_name' ,'Value', 'projection_x_coordinate'); % standard name
nc(ifld).Attribute(4) = struct('Name', 'comment' ,'Value', '1 km2 pixel centers');
ifld = ifld + 1;
nc(ifld).Name = 'y';
nc(ifld).Nctype = 'int';
nc(ifld).Dimension = {'y'};
nc(ifld).Attribute(1) = struct('Name', 'long_name' ,'Value', 'y-coordinate in Cartesian system');
nc(ifld).Attribute(2) = struct('Name', 'units' ,'Value', 'km');
nc(ifld).Attribute(3) = struct('Name', 'standard_name' ,'Value', 'projection_y_coordinate'); % standard name
nc(ifld).Attribute(4) = struct('Name', 'comment' ,'Value', '1 km2 pixel centers');
ifld = ifld + 1;
nc(ifld).Name = 'x_cor';
nc(ifld).Nctype = 'int';
nc(ifld).Dimension = {'x_cor'};
nc(ifld).Attribute(1) = struct('Name', 'long_name' ,'Value', 'x-coordinate of pixel corners in Cartesian system');
nc(ifld).Attribute(2) = struct('Name', 'units' ,'Value', 'km');
nc(ifld).Attribute(3) = struct('Name', 'standard_name' ,'Value', 'projection_x_coordinate'); % standard name
nc(ifld).Attribute(4) = struct('Name', 'comment' ,'Value', '1 km2 pixel corners');
ifld = ifld + 1;
nc(ifld).Name = 'y_cor';
nc(ifld).Nctype = 'int';
nc(ifld).Dimension = {'y_cor'};
nc(ifld).Attribute(1) = struct('Name', 'long_name' ,'Value', 'y-coordinate of pixel corners in Cartesian system');
nc(ifld).Attribute(2) = struct('Name', 'units' ,'Value', 'km');
nc(ifld).Attribute(3) = struct('Name', 'standard_name' ,'Value', 'projection_y_coordinate'); % standard name
nc(ifld).Attribute(4) = struct('Name', 'comment' ,'Value', '1 km2 pixel corners');
if OPT.ll
%% Longitude
% http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html#longitude-coordinate
ifld = ifld + 1;
nc(ifld).Name = 'longitude';
nc(ifld).Nctype = 'float';
nc(ifld).Dimension = {'x','y'};
nc(ifld).Attribute(1) = struct('Name', 'long_name' ,'Value', 'longitude');
nc(ifld).Attribute(2) = struct('Name', 'units' ,'Value', 'degrees_east');
nc(ifld).Attribute(3) = struct('Name', 'standard_name' ,'Value', 'longitude'); % standard name
nc(ifld).Attribute(4) = struct('Name', 'actual_range' ,'Value', [min(D.loncen(:)) max(D.loncen(:))]); %
%% Latitude
% http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html#latitude-coordinate
ifld = ifld + 1;
nc(ifld).Name = 'latitude';
nc(ifld).Nctype = 'float';
nc(ifld).Dimension = {'x','y'};
nc(ifld).Attribute(1) = struct('Name', 'long_name' ,'Value', 'latitude');
nc(ifld).Attribute(2) = struct('Name', 'units' ,'Value', 'degrees_north');
nc(ifld).Attribute(3) = struct('Name', 'standard_name' ,'Value', 'latitude'); % standard name
nc(ifld).Attribute(4) = struct('Name', 'actual_range' ,'Value', [min(D.latcen(:)) max(D.latcen(:))]); %
end % if OPT.ll
%% Time
% http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html#time-coordinate
% time is a dimension, so there are two options:
% * the variable name needs the same as the dimension
% http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html#id2984551
% * there needs to be an indirect mapping through the coordinates attribute
% http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html#id2984605
ifld = ifld + 1;
nc(ifld).Name = 'time';
nc(ifld).Nctype = 'double'; % float not sufficient as datenums are big: double
nc(ifld).Dimension = {'time'}; % {'locations','time'} % does not work in ncBrowse, nor in Quickplot (is indirect time mapping)
nc(ifld).Attribute(1) = struct('Name', 'long_name' ,'Value', 'time');
nc(ifld).Attribute(2) = struct('Name', 'units' ,'Value', ['days since ',datestr(OPT.refdatenum,'yyyy-mm-dd'),' 00:00:00 ',OPT.timezone]);
nc(ifld).Attribute(3) = struct('Name', 'standard_name' ,'Value', 'time');
nc(ifld).Attribute(4) = struct('Name', '_FillValue' ,'Value', OPT.fillvalue);
%% Parameters with standard names
% * http://cf-pcmdi.llnl.gov/documents/cf-standard-names/standard-name-table/current/
%% Define dimensions in this order:
% time,z,y,x
ifld = ifld + 1;
nc(ifld).Name = 'SST';
nc(ifld).Nctype = 'double';
nc(ifld).Dimension = {'x','y','time'};
nc(ifld).Attribute( 1) = struct('Name', 'long_name' ,'Value', 'sea surface temperature');
nc(ifld).Attribute(end+1) = struct('Name', 'units' ,'Value', 'degrees_Celcius');
nc(ifld).Attribute(end+1) = struct('Name', 'standard_name' ,'Value', 'sea_surface_skin_temperature'); % standard name
nc(ifld).Attribute(end+1) = struct('Name', '_FillValue' ,'Value', OPT.fillvalue);
nc(ifld).Attribute(end+1) = struct('Name', 'coordinates' ,'Value', 'latitude longitude');
nc(ifld).Attribute(end+1) = struct('Name', 'grid_mapping' ,'Value', 'polar_stereographic');
nc(ifld).Attribute(end+1) = struct('Name', 'actual_range' ,'Value', [D.data_min_value D.data_max_value]);
if OPT.pack
nc(ifld).Nctype = 'int'; %'byte'; %short
nc(ifld).Attribute(9) = struct('Name', 'valid_min' ,'Value', D.count_min_value);
nc(ifld).Attribute(10) = struct('Name', 'valid_max' ,'Value', D.count_max_value);
nc(ifld).Attribute(11) = struct('Name', 'scale_factor' ,'Value', D.gain);
nc(ifld).Attribute(12) = struct('Name', 'add_offset' ,'Value', D.offset);
end
%nc(ifld).Attribute(6) = struct('Name', 'cell_methods' ,'Value', 'point');x
%nc(ifld).Attribute(6) = struct('Name', 'cell_methods' ,'Value', 'point');y
%nc(ifld).Attribute(6) = struct('Name', 'cell_methods' ,'Value', 'point');time
%% 4 Create variables with attibutes
for ifld=1:length(nc)
nc_addvar(outputfile, nc(ifld));
end
%% 5 Fill variables
nc_varput(outputfile, 'x' , [1:D.nx]'-0.5);
nc_varput(outputfile, 'y' , [1:D.ny]'-0.5);
nc_varput(outputfile, 'x_cor' , [1:(D.nx+1)]');
nc_varput(outputfile, 'y_cor' , [1:(D.ny+1)]');
if OPT.ll
nc_varput(outputfile, 'longitude', D.loncen');
nc_varput(outputfile, 'latitude' , D.latcen');
end % if OPT.ll
nc_varput(outputfile, 'time' , D.datenum' - OPT.refdatenum);
if OPT.pack
nc_datput(outputfile, 'SST' , int16(flipud(D.count)')); % do not use nc_varput with scale_factor and add_offset
% uint8 and int8 both becomes nc_byte,
% * upper D.count levels ( > 512) are outside n_byte range
% * unit 8 is read as int8, so [-512 512] instad of [0 1024]
% * altarnative int16 contains too much space (twice).
else
nc_varput(outputfile, 'SST' , D.data');
end
%% 6 Check
if OPT.dump
nc_dump(outputfile);
end
%% Pause
if OPT.pause
pausedisp
end
end %for ifile=1:length(OPT.files)