WSL/SLF GitLab Repository

NetCDFIO.cc 49.4 KB
Newer Older
1
/***********************************************************************************/
2
/*  Copyright 2014 WSL Institute for Snow and Avalanche Research    SLF-DAVOS      */
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
/***********************************************************************************/
/* This file is part of MeteoIO.
    MeteoIO is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    MeteoIO is distributed 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.  See the
    GNU Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public License
    along with MeteoIO.  If not, see <http://www.gnu.org/licenses/>.
*/
#include "NetCDFIO.h"
19
#include <meteoio/ResamplingAlgorithms2D.h>
20
#include <meteoio/meteoStats/libinterpol1D.h>
21
#include <meteoio/Timer.h>
22
#include <meteoio/MathOptim.h>
23
24
25
26
27
#include <meteoio/plugins/libncpp.h>

#include <cmath>
#include <cstdio>
#include <algorithm>
28
29

using namespace std;
30
using namespace ncpp; // wrappers for libnetcdf
31
32
33
34
35

namespace mio {
/**
 * @page netcdf NetCDF
 * @section netcdf_format Format
36
37
38
39
40
 * In order to promote creation, access and sharing of scientific data, the NetCDF format has been
 * created as a machine-independent format. NetCDF (network Common Data Form) is therefore an interface
 * for array-oriented data access and a library that provides an implementation of the interface. The
 * <A HREF="http://www.unidata.ucar.edu/downloads/netcdf/index.jsp">NetCDF software</A> was developed
 * at the <A HREF="http://www.unidata.ucar.edu/">Unidata Program Center</A> in Boulder, Colorado.
41
42
 * In order to graphicaly explore the content and structure of NetCDF files, you can use the
 * <A REF="http://www.epic.noaa.gov/java/ncBrowse/">ncBrowse></A> java software.
43
44
45
46
47
48
49
50
 *
 * The <A HREF="http://cfconventions.org/1.6.html">conventions</A> for climate and forecast (CF) metadata
 * are designed to promote the processing and sharing of netCDF files. The conventions define metadata
 * that provide a definitive description of what the data represents, and the spatial and temporal properties of the data.
 * This plugin follows such conventions as well as the naming extensions defined by the
 * <A HREF="http://www.cnrm.meteo.fr/">CNRM</A>.
 *
 * *Put here the more informations about the standard format that is implemented*
51
52
53
54
55
56
57
58
 *
 * @section netcdf_units Units
 *
 *
 * @section netcdf_keywords Keywords
 * This plugin uses the following keywords:
 * - COORDSYS: coordinate system (see Coords); [Input] and [Output] section
 * - COORDPARAM: extra coordinates parameters (see Coords); [Input] and [Output] section
59
60
 * - DEMFILE: The filename of the file containing the DEM; [Input] section
 * - DEMVAR: The variable name of the DEM within the DEMFILE; [Input] section
61
 * - METEOFILE: the NetCDF file which shall be used for the meteo parameter input/output; [Input] and [Output] section
62
 * - GRID2DFILE: the NetCDF file which shall be used for gridded input/output; [Input] and [Output] section
63
64
 * - STRICTFORMAT: Whether the NetCDF file should be strictly compliant with the CNRM standard; Parameters not present
 *                 in the specification will be omitted; [Input] and [Output] section
65
 *
66
 * @section netcdf_example Example use
67
68
69
70
71
72
73
 * @code
 * [Input]
 * DEM     = NETCDF
 * DEMFILE = ./input/Aster_tile.nc
 * DEMVAR  = z
 * @endcode
 *
74
 * @section netcdf_compilation Compilation
75
76
 * In order to compile this plugin, you need libnetcdf (for C). For Linux, please select both the libraries and
 * their development files in your package manager.
77
78
 */

79
const double NetCDFIO::plugin_nodata = -9999999.; //CNRM-GAME nodata value
80
const double NetCDFIO::epsilon = 1.0e-10; //when comparing timestamps
81

82
83
84
const std::string NetCDFIO::cf_time = "time";
const std::string NetCDFIO::cf_units = "units";
const std::string NetCDFIO::cf_days = "days since ";
85
const std::string NetCDFIO::cf_hours = "hours since ";
86
const std::string NetCDFIO::cf_seconds = "seconds since ";
87
88
89
90
91
92
const std::string NetCDFIO::cf_latitude = "lat";
const std::string NetCDFIO::cf_longitude = "lon";
const std::string NetCDFIO::cf_altitude = "z";
const std::string NetCDFIO::cf_ta = "temperature";
const std::string NetCDFIO::cf_rh = "humidity";
const std::string NetCDFIO::cf_p = "pressure";
93

94
95
96
const std::string NetCDFIO::cnrm_points = "Number_of_points";
const std::string NetCDFIO::cnrm_latitude = "LAT";
const std::string NetCDFIO::cnrm_longitude = "LON";
97
98
99
const std::string NetCDFIO::cnrm_altitude = "ZS";
const std::string NetCDFIO::cnrm_aspect = "aspect";
const std::string NetCDFIO::cnrm_slope = "slope";
100
101
102
103
104
105
106
107
108
109
110
111
112
113
const std::string NetCDFIO::cnrm_ta = "Tair";
const std::string NetCDFIO::cnrm_rh = "HUMREL";
const std::string NetCDFIO::cnrm_vw = "Wind";
const std::string NetCDFIO::cnrm_dw = "Wind_DIR";
const std::string NetCDFIO::cnrm_qair = "Qair";
const std::string NetCDFIO::cnrm_co2air = "CO2air";
const std::string NetCDFIO::cnrm_theorsw = "theorSW";
const std::string NetCDFIO::cnrm_neb = "NEB";
const std::string NetCDFIO::cnrm_hnw = "Rainf";
const std::string NetCDFIO::cnrm_snowf = "Snowf";
const std::string NetCDFIO::cnrm_swr_direct = "DIR_SWdown";
const std::string NetCDFIO::cnrm_swr_diffuse = "SCA_SWdown";
const std::string NetCDFIO::cnrm_p = "PSurf";
const std::string NetCDFIO::cnrm_ilwr = "LWdown";
Thomas Egger's avatar
Thomas Egger committed
114
const std::string NetCDFIO::cnrm_timestep = "FRC_TIME_STP";
115

116
117
118
119
120
121
122
123
124
const std::string NetCDFIO::ecmwf_ta = "t2m";
const std::string NetCDFIO::ecmwf_p = "sp";
const std::string NetCDFIO::ecmwf_iswr = "ssrd";
const std::string NetCDFIO::ecmwf_ilwr = "strd";
const std::string NetCDFIO::ecmwf_hnw = "tp";
const std::string NetCDFIO::ecmwf_td = "d2m";
const std::string NetCDFIO::ecmwf_u10 = "u10";
const std::string NetCDFIO::ecmwf_v10 = "v10";

125
std::map<std::string, size_t> NetCDFIO::paramname;
126
std::map<std::string, std::string> NetCDFIO::map_name;
127
128
129
130
131
const bool NetCDFIO::__init = NetCDFIO::initStaticData();

bool NetCDFIO::initStaticData()
{
	//Associate unsigned int value and a string representation of a meteo parameter
132
	paramname[cnrm_ta] = MeteoData::TA;
133
134
135
136
	paramname[cnrm_qair] = IOUtils::npos; // not a standard MeteoIO parameter
	paramname[cnrm_co2air] = IOUtils::npos; // not a standard MeteoIO parameter
	paramname[cnrm_neb] = IOUtils::npos; // not a standard MeteoIO parameter
	paramname[cnrm_theorsw] = IOUtils::npos; // not a standard MeteoIO parameter
137
138
139
140
141
142
143
144
145
	paramname[cnrm_rh] = MeteoData::RH;
	paramname[cnrm_vw] = MeteoData::VW;
	paramname[cnrm_dw] = MeteoData::DW;
	paramname[cnrm_hnw] = IOUtils::npos;
	paramname[cnrm_snowf] = IOUtils::npos;
	paramname[cnrm_swr_direct] = IOUtils::npos;
	paramname[cnrm_swr_diffuse] = IOUtils::npos;
	paramname[cnrm_p] = MeteoData::P;
	paramname[cnrm_ilwr] = MeteoData::ILWR;
146

147
148
149
150
151
152
	map_name["TA"] = cnrm_ta;
	map_name["RH"] = cnrm_rh;
	map_name["ILWR"] = cnrm_ilwr;
	map_name["P"] = cnrm_p;
	map_name["VW"] = cnrm_vw;
	map_name["DW"] = cnrm_dw;
Thomas Egger's avatar
Thomas Egger committed
153
	map_name["ISWR"] = cnrm_swr_direct;
154
	map_name["HNW"] = cnrm_hnw;
155
156
157
158
	map_name[cnrm_co2air] = cnrm_co2air;
	map_name[cnrm_qair] = cnrm_qair;
	map_name[cnrm_theorsw] = cnrm_theorsw;
	map_name[cnrm_neb] = cnrm_neb;
159

160
161
162
	return true;
}

163
NetCDFIO::NetCDFIO(const std::string& configfile) : cfg(configfile), coordin(), coordinparam(), coordout(), coordoutparam(),
164
                                                    in_dflt_TZ(0.), out_dflt_TZ(0.), in_strict(false), out_strict(false), vecMetaData()
165
166
{
	IOUtils::getProjectionParameters(cfg, coordin, coordinparam, coordout, coordoutparam);
167
	parseInputOutputSection();
168
169
}

170
NetCDFIO::NetCDFIO(const Config& cfgreader) : cfg(cfgreader), coordin(), coordinparam(), coordout(), coordoutparam(),
171
                                              in_dflt_TZ(0.), out_dflt_TZ(0.), in_strict(false), out_strict(false), vecMetaData()
172
173
{
	IOUtils::getProjectionParameters(cfg, coordin, coordinparam, coordout, coordoutparam);
174
	parseInputOutputSection();
175
176
}

177
NetCDFIO::~NetCDFIO() throw() {}
178

179
180
181
182
183
184
void NetCDFIO::parseInputOutputSection()
{
	//default timezones
	in_dflt_TZ = out_dflt_TZ = IOUtils::nodata;
	cfg.getValue("TIME_ZONE", "Input", in_dflt_TZ, IOUtils::nothrow);
	cfg.getValue("TIME_ZONE", "Output", out_dflt_TZ, IOUtils::nothrow);
185
186
187

	cfg.getValue("STRICTFORMAT", "Input", in_strict, IOUtils::nothrow);
	cfg.getValue("STRICTFORMAT", "Output", out_strict, IOUtils::nothrow);
188
189
}

190
void NetCDFIO::read2DGrid(Grid2DObject& grid_out, const std::string& arguments)
191
{
192
193
194
195
196
197
198
199
	vector<string> vec_argument;
	IOUtils::readLineToVec(arguments, vec_argument, ':');

	if (vec_argument.size() == 2) {
		read2DGrid_internal(grid_out, vec_argument[0], vec_argument[1]);
	} else {
		throw InvalidArgumentException("The format for the arguments to NetCDFIO::read2DGrid is filename:varname", AT);
	}
200
201
}

202
void NetCDFIO::read2DGrid(Grid2DObject& grid_out, const MeteoGrids::Parameters& parameter, const Date& date)
203
{
204
	string filename;
205
206
	cfg.getValue("GRID2DFILE", "Input", filename);

207
	const string varname = get_varname(parameter);
208

209
	read2DGrid_internal(grid_out, filename, varname, date);
210
211
}

212
void NetCDFIO::read2DGrid_internal(Grid2DObject& grid_out, const std::string& filename, const std::string& varname, const Date& date)
213
{
214
	const bool is_record = (date != Date());
215
216
	size_t lat_index = 0, lon_index = 1;

217
	int ncid, varid;
218
	vector<int> dimid, dim_varid;
219
220
221
222
223
224
225
	vector<string> dimname;
	vector<size_t> dimlen;

	open_file(filename, NC_NOWRITE, ncid);
	get_variable(ncid, varname, varid);
	get_dimension(ncid, varname, varid, dimid, dim_varid, dimname, dimlen);

226
	if (is_record) { // In case we're reading a record the first index is always the record index
227
228
229
230
231
		lat_index = 1;
		lon_index = 2;

		if (dimid.size()!=3 || dimlen[0]<1 || dimlen[lat_index]<2 || dimlen[lon_index]<2)
			throw IOException("Variable '" + varname + "' may only have three dimensions, all have to at least have length 1", AT);
232
233
234
235
236
237
	} else if (dimid.size()==3 && dimlen[0]==1) { //in case the variable is associated with a 1 element time dimension
		lat_index = 1;
		lon_index = 2;

		if (dimlen[lat_index]<2 || dimlen[lon_index]<2)
			throw IOException("All dimensions for variable '" + varname + "' have to at least have length 1", AT);
238
	} else if (dimid.size()!=2 || dimlen[lat_index]<2 || dimlen[lon_index]<2) {
239
		throw IOException("Variable '" + varname + "' may only have two dimensions and both have to have length >1", AT);
240
241
242
243
244
	}

	double *lat = new double[dimlen[lat_index]];
	double *lon = new double[dimlen[lon_index]];
	double *grid = new double[dimlen[lat_index]*dimlen[lon_index]];
245

246
247
	read_data(ncid, dimname[lat_index], dim_varid[lat_index], lat);
	read_data(ncid, dimname[lon_index], dim_varid[lon_index], lon);
248

249
	if (is_record) {
250
		const size_t pos = find_record(ncid, NetCDFIO::cf_time, dimid[0], date.getModifiedJulianDate());
251
252
		if (pos == IOUtils::npos)
			throw IOException("No record for date " + date.toString(Date::ISO), AT);
253

254
255
256
257
		read_data(ncid, varname, varid, pos, dimlen[lat_index], dimlen[lon_index], grid);
	} else {
		read_data(ncid, varname, varid, grid);
	}
258

259
260
261
	double missing_value=plugin_nodata;
	if (ncpp::check_attribute(ncid, varid, "missing_value"))
		ncpp::get_attribute(ncid, varname, varid, "missing_value", missing_value);
262

263
264
265
266
267
268
269
270
271
272
273
274
275
	copy_grid(dimlen[lat_index], dimlen[lon_index], lat, lon, grid, missing_value, grid_out);

	//handle data packing if necessary
	if (ncpp::check_attribute(ncid, varid, "scale_factor")) {
		double scale_factor=1.;
		ncpp::get_attribute(ncid, varname, varid, "scale_factor", scale_factor);
		grid_out.grid2D *= scale_factor;
	}
	if (ncpp::check_attribute(ncid, varid, "add_offset")) {
		double add_offset=0.;
		ncpp::get_attribute(ncid, varname, varid, "add_offset", add_offset);
		grid_out.grid2D += add_offset;
	}
276

277
	close_file(filename, ncid);
278
	delete[] lat; delete[] lon; delete[] grid;
279
280
}

281
void NetCDFIO::copy_grid(const size_t& latlen, const size_t& lonlen, const double * const lat, const double * const lon,
282
                         const double * const grid, const double& nodata, Grid2DObject& grid_out)
283
{
284
	double resampling_factor_x = IOUtils::nodata, resampling_factor_y=IOUtils::nodata;
285
	const double cellsize = calculate_cellsize(latlen, lonlen, lat, lon, resampling_factor_x, resampling_factor_y);
286

287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
	const double cntr_lat = .5*(lat[0]+lat[latlen-1]);
	const double cntr_lon = .5*(lon[0]+lon[lonlen-1]);

	//computing lower left corner by using the center point as reference
	Coords cntr(coordin, coordinparam);
	cntr.setLatLon(cntr_lat, cntr_lon, IOUtils::nodata);
	cntr.moveByXY(-.5*(double)(lonlen-1)*cellsize, -.5*(double)(latlen-1)*cellsize);

	grid_out.set(lonlen, latlen, cellsize, cntr);

	//Handle the case of llcorner/urcorner swapped
	if (lat[0]<=lat[latlen-1]) {
		for (size_t kk=0; kk < latlen; kk++) {
			const size_t row = kk*lonlen;
			if (lon[0]<=lon[lonlen-1]) {
				for (size_t ll=0; ll < lonlen; ll++)
					grid_out(ll, kk) = IOUtils::standardizeNodata(grid[row + ll], nodata);
			} else {
				for (size_t ll=0; ll < lonlen; ll++)
					grid_out(ll, kk) = IOUtils::standardizeNodata(grid[row + (lonlen -1) - ll], nodata);
			}
		}
	} else {
		for (size_t kk=0; kk < latlen; kk++) {
			const size_t row = ((latlen-1) - kk)*lonlen;
			if (lon[0]<=lon[lonlen-1]) {
				for (size_t ll=0; ll < lonlen; ll++)
					grid_out(ll, kk) = IOUtils::standardizeNodata(grid[row + ll], nodata);
			} else {
				for (size_t ll=0; ll < lonlen; ll++)
					grid_out(ll, kk) = IOUtils::standardizeNodata(grid[row + (lonlen -1) - ll], nodata);
			}
319
320
321
		}
	}

322
323
	if (resampling_factor_x != IOUtils::nodata) {
		grid_out.grid2D = ResamplingAlgorithms2D::BilinearResampling(grid_out.grid2D, resampling_factor_x, resampling_factor_y);
324
325
326
		grid_out.ncols = grid_out.grid2D.getNx();
		grid_out.nrows = grid_out.grid2D.getNy();
	}
327
328
}

329
330
331
/* The Grid2DObject holds data and meta data for quadratic cells. However the NetCDF file
 * stores the grid as discrete latitude and longitude values. It is necessary to calculate
 * the distance between the edges of the grid and determine the cellsize. This cellsize may
332
333
 * be different for X and Y directions. We then choose one cellsize for our grid and
 * determine a factor that will be used for resampling the grid to likewise consist of
334
335
 * quadratic cells.
 */
336
337
double NetCDFIO::calculate_cellsize(const size_t& latlen, const size_t& lonlen, const double * const lat, const double * const lon,
                                    double& factor_x, double& factor_y)
338
{
339
340
341
342
343
344
	const double cntr_lat = .5*(lat[0]+lat[latlen-1]);
	const double cntr_lon = .5*(lon[0]+lon[lonlen-1]);
	double alpha;

	const double distanceX = Coords::VincentyDistance(cntr_lat, lon[0], cntr_lat, lon[lonlen-1], alpha);
	const double distanceY = Coords::VincentyDistance(lat[0], cntr_lon, lat[latlen-1], cntr_lon, alpha);
345

346
	// lonlen, latlen are decremented by 1; n linearly connected points have (n-1) connections
347
348
	const double cellsize_x = distanceX / (lonlen-1);
	const double cellsize_y = distanceY / (latlen-1);
349

350
351
	// round to 1cm precision for numerical stability
	const double cellsize = Optim::round( std::min(cellsize_x, cellsize_y)*100. ) / 100.;
352

353
354
	if (cellsize_x == cellsize_y) {
		return cellsize_x;
355
	} else {
356
357
358
359
		factor_x =  cellsize_x / cellsize;
		factor_y =  cellsize_y / cellsize;

		return cellsize;
360
	}
361
362
}

363
364
void NetCDFIO::readDEM(DEMObject& dem_out)
{
365
	string filename, varname;
366

367
368
	cfg.getValue("DEMFILE", "Input", filename);
	cfg.getValue("DEMVAR", "Input", varname);
369

370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
	read2DGrid_internal(dem_out, filename, varname);
}

void NetCDFIO::readLanduse(Grid2DObject& /*landuse_out*/)
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

void NetCDFIO::readAssimilationData(const Date& /*date_in*/, Grid2DObject& /*da_out*/)
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

void NetCDFIO::readStationData(const Date&, std::vector<StationData>& vecStation)
{
387
	if (!vecMetaData.empty()) { // We already have meta data
388
389
390
391
		vecStation = vecMetaData;
		return;
	}

392
	string filename;
393
394
395
396
397
398
399
400
401
402
403
404
	cfg.getValue("METEOFILE", "Input", filename);

	int ncid;
	open_file(filename, NC_NOWRITE, ncid);
	readMetaData(ncid, vecMetaData);
	close_file(filename, ncid);

	vecStation = vecMetaData;
}

void NetCDFIO::readMetaData(const int& ncid, std::vector<StationData>& vecStation)
{
405
406
407
	vecStation.clear();

	int dimid;
408
	size_t dimlen;
409
	map<string, int> map_vid;
410

411
	get_dimension(ncid, cnrm_points, dimid, dimlen);
412
413
414
	if (dimlen == 0) return; // There are no stations

	get_meta_data_ids(ncid, map_vid);
415
416
417
418

	double *alt = new double[dimlen];
	double *lat = new double[dimlen];
	double *lon = new double[dimlen];
419
420
	double *aspect = new double[dimlen];
	double *slope = new double[dimlen];
421

422
423
424
425
426
	read_data(ncid, cnrm_altitude, map_vid[cnrm_altitude], alt);
	read_data(ncid, cnrm_latitude, map_vid[cnrm_latitude], lat);
	read_data(ncid, cnrm_longitude, map_vid[cnrm_longitude], lon);
	read_data(ncid, cnrm_aspect, map_vid[cnrm_aspect], aspect);
	read_data(ncid, cnrm_slope, map_vid[cnrm_slope], slope);
427
428
429
430
431
432

	//Parse to StationData objects
	Coords location(coordin, coordinparam);
	ostringstream ss;
	for (size_t ii=0; ii<dimlen; ii++) {
		location.setLatLon(lat[ii], lon[ii], alt[ii]);
433

434
		ss << (ii+1);
435
		const string id( ss.str() );
436
437
438
		ss.str("");

		ss << "Station " << (ii +1);
439
		const string name( ss.str() );
440
		ss.str("");
441

442
		StationData tmp(location, id, name);
443
		const double aspect_bearing = (aspect[ii] < 0) ? 0 : aspect[ii]; // aspect allowed to be -1 in CNRM format...
444
		tmp.setSlope(slope[ii], aspect_bearing);
445
446
447
		vecStation.push_back(tmp);
	}

448
	delete[] alt; delete[] lat; delete[] lon; delete[] aspect; delete[] slope;
449
450
}

451
452
453
454
455
456
457
458
void NetCDFIO::get_meta_data_ids(const int& ncid, std::map<std::string, int>& map_vid)
{
	const string names[] = {cnrm_altitude, cnrm_latitude, cnrm_longitude, cnrm_aspect, cnrm_slope};
	vector<string> varname(names, names + sizeof(names) / sizeof(names[0]));

	vector<string> dimensions;
	dimensions.push_back(cnrm_points); // All variables have to have the dimension cnrm_points

459
	for (vector<string>::const_iterator it = varname.begin(); it != varname.end(); ++it) {
460
461
462
463
		int varid;
		const string& name = *it;

		get_variable(ncid, name, varid);
464
465
		if (!check_dimensions(ncid, name, varid, dimensions))
			throw IOException("Variable '" + name  + "' fails dimension check", AT);
466
467
468
469
470

		map_vid[name] = varid;
	}
}

471
void NetCDFIO::readMeteoData(const Date& dateStart, const Date& dateEnd, std::vector< std::vector<MeteoData> >& vecMeteo, const size_t&)
472
{
473
474
	vecMeteo.clear();

475
	string filename;
476
477
478
479
480
481
482
483
484
485
486
	cfg.getValue("METEOFILE", "Input", filename);

	int ncid;
	open_file(filename, NC_NOWRITE, ncid);

	if (vecMetaData.empty()) readMetaData(ncid, vecMetaData);

	if (!vecMetaData.empty()) { //at least one station exists
		size_t index_start, index_end;
		vector<Date> vec_date;
		get_indices(ncid, dateStart, dateEnd, index_start, index_end, vec_date); //get indices for dateStart and dateEnd
487
488

		MeteoData meteo_data; //the template MeteoData object
489
		if ((index_start != IOUtils::npos) && (index_end != IOUtils::npos)) {
490
			map<string, size_t> map_parameters;
491
			get_parameters(ncid, map_parameters, meteo_data); //get a list of parameters present an render the template
492

493
			readData(ncid, index_start, vec_date, map_parameters, meteo_data, vecMeteo);
494
495
496
497
498
499
		}
	}

	close_file(filename, ncid);
}

500
void NetCDFIO::readData(const int& ncid, const size_t& index_start, const std::vector<Date>& vec_date,
501
                        const std::map<std::string, size_t>& map_parameters, const MeteoData& meteo_data, std::vector< std::vector<MeteoData> >& vecMeteo)
502
{
503
504
	const size_t number_of_stations = vecMetaData.size();
	const size_t number_of_records = vec_date.size();
505

506
	// Allocate all the MeteoData objects based on the template meteo_data
507
	vector<MeteoData> tmp_vec(number_of_records, meteo_data);
508
509
510
511
512
513
514
	for (size_t jj=0; jj<number_of_records; jj++) tmp_vec[jj].date = vec_date[jj]; //set correct date for every record

	for (size_t ii=0; ii<number_of_stations; ii++) {
		for (size_t jj=0; jj<number_of_records; jj++) tmp_vec[jj].meta = vecMetaData[ii]; //adapt meta data
		vecMeteo.push_back(tmp_vec);
	}

515
	// Allocate enough linear space for each parameter and read the data from NetCDF
516
	map<string, double*> map_data;
517
	for (map<string, size_t>::const_iterator it = map_parameters.begin(); it != map_parameters.end(); ++it) {
518
		double* data = new double[number_of_stations*number_of_records];
519
		const string& varname = it->first;
520

521
		map_data[varname] = data;
522

523
		int varid;
524
525
526
527
		get_variable(ncid, varname, varid);
		read_data_2D(ncid, varname, varid, index_start, number_of_records, number_of_stations, data);
	}

528
	copy_data(ncid, map_parameters, map_data, number_of_stations, number_of_records, vecMeteo);
529

530
	for (map<string, double*>::const_iterator it = map_data.begin(); it != map_data.end(); ++it) {
531
532
533
534
		delete[] it->second;
	}
}

535
536
537
538
539
540
// The copying of data into vecMeteo is a process consisting of:
// 1. A check what the relation between MeteoIO parameters and CNRM parameters present is, check map_parameters
// 2. If there is no direct association between the parameters present and the meteo_data parameters we might
//    have to deal with the parameter in a more complex way: e.g., HNW or SWR measurements
// 3. Once we know how to deal with the parameter we loop through all stations and all parameters and copy them
//    into the appropriate places. All unit conversion have been accomplished at that point.
541
void NetCDFIO::copy_data(const int& ncid, const std::map<std::string, size_t>& map_parameters, const std::map<std::string, double*> map_data,
542
                         const size_t& number_of_stations, const size_t& number_of_records, std::vector< std::vector<MeteoData> >& vecMeteo)
543
{
544
	for (map<string, double*>::const_iterator it = map_data.begin(); it != map_data.end(); ++it) {
545
546
547
		const string& varname = it->first;

		//find correct handling for each parameter
548
		bool simple_copy = false, mutiply_copy = false, hnw_measurement = false, sw_measurement = false;
549
		double multiplier = IOUtils::nodata;
550
		const size_t param = map_parameters.find(varname)->second; //must exist, at this point we know it does
551
552

		if (param == IOUtils::npos) {
553
554
			if ((varname == cnrm_snowf) || (varname == cnrm_hnw)) {
				int varid;
Thomas Egger's avatar
Thomas Egger committed
555
556
				get_variable(ncid, cnrm_timestep, varid);
				read_value(ncid, cnrm_timestep, varid, multiplier);
557

Thomas Egger's avatar
Thomas Egger committed
558
				if (multiplier <= 0) throw InvalidArgumentException("The variable '" + cnrm_timestep + "' is invalid", AT);
559
560
561
562

				hnw_measurement = true;
			} else if ((varname == cnrm_swr_diffuse) || (varname == cnrm_swr_direct)) {
				sw_measurement = true;
563
564
			} else {
				throw IOException("Don't know how to deal with parameter " + varname, AT);
565
566
			}
		} else {
567
			if (varname == cnrm_rh) {
568
569
570
571
572
573
574
				mutiply_copy = true;
				multiplier = 0.01;
			} else {
				simple_copy = true;
			}
		}

575
		// Loop through all times and all stations
576
		for (size_t jj=0; jj<number_of_records; jj++) {
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
			for (size_t ii=0; ii<number_of_stations; ii++) {
				double& value = (it->second)[jj*number_of_stations + ii];
				bool nodata = false;

				if (value == plugin_nodata) {
					nodata = true;
					value = IOUtils::nodata;
				}

				if (simple_copy) {
					vecMeteo[ii][jj](param) = value;
				} else if (mutiply_copy) {
					if (nodata) {
						vecMeteo[ii][jj](param) = value;
					} else {
						vecMeteo[ii][jj](param) = value * multiplier;
					}
594
				} else if (hnw_measurement) {
595
					if (!nodata) {
596
						double& hnw = vecMeteo[ii][jj](MeteoData::HNW);
597
598
599
						if (hnw == IOUtils::nodata) hnw = 0.0;
						hnw += value * multiplier;
					}
600
601
602
603
				} else if (sw_measurement) {
					if (!nodata) {
						double& iswr = vecMeteo[ii][jj](MeteoData::ISWR);
						if (iswr == IOUtils::nodata) iswr = 0.0;
604
						iswr += value;
605
					}
606
				}
607
608
609
610
611
			}
		}
	}
}

612
613
614
// Go through all variables present in the NetCDF dataset that have the correct dimensions. A map called
// map_parameters will associate all parameters present with MeteoData parameters or IOUtils::npos). If
// the CNRM parameter does not have a corresponding parameter in the meteo_data object we can add a new
615
616
// parameter (e.g. cnrm_theorsw) or if the situation is more complex (e.g. rainfall is measured with two
// parameters) we deal with the situation in copy_data().
617
void NetCDFIO::get_parameters(const int& ncid, std::map<std::string, size_t>& map_parameters, MeteoData& meteo_data)
618
{
619
620
	vector<string> dimensions;
	dimensions.push_back(cf_time);
621
622
623
624
625
626
627
	dimensions.push_back(cnrm_points);

	vector<string> present_parameters;
	get_variables(ncid, dimensions, present_parameters);

	for (vector<string>::const_iterator it = present_parameters.begin(); it != present_parameters.end(); ++it) {
		const string& name = *it;
628
		//cout << "Found parameter: " << name << endl;
629

630
		// Check if parameter exists in paramname, which holds strict CNRM parameters
631
		const map<string, size_t>::const_iterator strict_it = paramname.find(name);
632
633
		if (strict_it != paramname.end()) { // parameter is a part of the CNRM specification
			size_t index = strict_it->second;
634
635
636
637
638

			if ((name == cnrm_theorsw) || (name == cnrm_qair) || (name == cnrm_co2air) || (name == cnrm_neb)) {
			 	index = meteo_data.addParameter(name);
			}

639
640
641
			map_parameters[name] = index;
		} else if (!in_strict) { // parameter will be read anyway
			size_t index = IOUtils::npos;
642

643
644
645
646
647
648
649
			if (meteo_data.param_exists(name)) {
				index = meteo_data.getParameterIndex(name);
			} else {
			 	index = meteo_data.addParameter(name);
			}

			map_parameters[name] = index;
650
651
652
653
		}
	}
}

654
655
656
657
// The CNRM format stores timestamps as doubles (either seconds or days counted from a start date)
// This method takes the dateStart and dateEnd requested and looks for the corresponding indices
// of the time variable indexStart and indexEnd.
// Furthermore the timestamps are converted to mio::Date objects and stored in vecDate
658
659
void NetCDFIO::get_indices(const int& ncid, const Date& dateStart, const Date& dateEnd, size_t& indexStart, size_t& indexEnd, std::vector<Date>& vecDate)
{
660
661
	indexStart = indexEnd = IOUtils::npos;

662
663
664
665
666
	int varid, dimid;
	size_t dimlen;
	get_dimension(ncid, NetCDFIO::cf_time, dimid, dimlen);
	get_variable(ncid, NetCDFIO::cf_time, varid);

667
	// Get the units attribute and calculate the offset date
668
669
670
	string units_str;
	NetCDFIO::TimeUnit unit_type;
	Date offset;
671
672
	get_attribute(ncid, NetCDFIO::cf_time, varid, cf_units, units_str);
	calculate_offset(units_str, unit_type, offset);
673
674
675
676

	double *time = new double[dimlen];
	read_data(ncid, NetCDFIO::cf_time, varid, time);

677
	// Firstly, check whether search makes any sense, that is dateStart and dateEnd overlap with the times present
678
679
680
	bool search = true;
	if (dimlen > 0) {
		Date time_start(offset), time_end(offset);
681

682
683
684
685
686
687
688
		double start = time[0];
		double end = time[dimlen-1];

		if (unit_type == seconds) {
			start /= 86400;
			end   /= 86400;
		}
689
690
691
692
		if (unit_type == hours) {
			start /= 24;
			end   /= 24;
		}
693
694
		time_start += Date(start, 0.0);
		time_end += Date(end, 0.0);
695
696
697
698
699

		if (time_start > dateEnd) search = false;
		if (time_end < dateStart) search = false;
	}

700
701
	// If search is feasible then loop through the existent timestamps and find the relevant indices
	bool start_found = false;
702
703
704
705
706
	if (search) {
		for (size_t ii=0; ii<dimlen; ii++) {
			if (unit_type == seconds) {
				time[ii] /= 86400;
			}
707
708
709
			if (unit_type == hours) {
				time[ii] /= 24;
			}
710

711
			const Date tmp_date = offset + Date(time[ii], 0.0);
712

713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
			if (!start_found && (dateStart <= tmp_date && tmp_date <= dateEnd)) {
				start_found = true;
				indexStart = ii;
			} else if (start_found && (tmp_date > dateEnd)) {
				indexEnd = ii-1;
				break;
			}

			if (start_found) vecDate.push_back(tmp_date);
		}

		if (start_found && (indexEnd == IOUtils::npos)) {
			indexEnd = dimlen-1;
		}
	}
728

729
730
731
	delete[] time;
}

732
733
// The CNRM timestamps have an offset that is saved in the units attribute of
// the time variable - this method retrieves that offset
734
735
736
void NetCDFIO::calculate_offset(const std::string& units, NetCDFIO::TimeUnit& time_unit, Date& offset)
{
	string tmp(units);
737
	const size_t found_sec = units.find(NetCDFIO::cf_seconds);
738
	const size_t found_hour = units.find(NetCDFIO::cf_hours);
739
	const size_t found_day = units.find(NetCDFIO::cf_days);
740
741
742
743

	if (found_sec != string::npos) {
		time_unit = seconds;
		tmp = tmp.substr(found_sec + NetCDFIO::cf_seconds.size());
744
745
746
	} else if (found_hour != string::npos) {
		time_unit = hours;
		tmp = tmp.substr(found_hour + NetCDFIO::cf_hours.size());
747
748
	} else if (found_day != string::npos) {
		time_unit = days;
749
		tmp = tmp.substr(found_day + NetCDFIO::cf_days.size());
750
	} else {
751
		throw InvalidFormatException("Variable '"+NetCDFIO::cf_time+"' has no valid attribute '" + cf_units + "'" , AT);
752
	}
753

754
	const bool success = IOUtils::convertString(offset, tmp, in_dflt_TZ);
755
	if (!success) throw InvalidFormatException("Cannot parse time: " + tmp, AT);
756
757
}

758
void NetCDFIO::writeMeteoData(const std::vector< std::vector<MeteoData> >& vecMeteo, const std::string&)
759
{
760
	const size_t number_of_stations = vecMeteo.size();
761
762
	if (number_of_stations == 0) return; //Nothing to write

763
	const size_t number_of_records = vecMeteo[0].size();
764

765
	string filename;
766
767
768
769
	cfg.getValue("METEOFILE", "Output", filename);

	int ncid, did_time, vid_time, did_points;
	bool create_time = false, create_points = false, create_locations = false, create_variables = false;
770

771
	const bool exists = IOUtils::fileExists(filename);
772
	if (exists) remove(filename.c_str()); // NOTE: file is deleted if it exists
773
774

	double* dates;
775
776
777
	map<string, double*> map_data; // holds a pointer for every C array to be written
	map_data[cnrm_latitude] = new double[number_of_stations];
	map_data[cnrm_longitude] = new double[number_of_stations];
Thomas Egger's avatar
Thomas Egger committed
778
779
780
	map_data[cnrm_altitude] = new double[number_of_stations];
	map_data[cnrm_aspect] = new double[number_of_stations];
	map_data[cnrm_slope] = new double[number_of_stations];
781
782
783
784

	map<string, int> varid;
	map<size_t, string> map_param_name;

Thomas Egger's avatar
Thomas Egger committed
785
	get_parameters(vecMeteo, map_param_name, map_data, dates);
786

787
788
	create_file(filename, NC_CLASSIC_MODEL, ncid);
	create_time = create_points = create_locations = create_variables = true;
789
790
791

	if (create_time) create_time_dimension(ncid, did_time, vid_time);
	if (create_points) add_dimension(ncid, cnrm_points, number_of_stations, did_points);
Thomas Egger's avatar
Thomas Egger committed
792
793
	if (create_locations) create_meta_data(ncid, did_points, map_data, varid);
	if (create_variables) create_parameters(ncid, did_time, did_points, number_of_records, number_of_stations, map_param_name, map_data, varid);
794
795
796

	end_definitions(filename, ncid);

Thomas Egger's avatar
Thomas Egger committed
797
	copy_data(number_of_stations, number_of_records, vecMeteo, map_param_name, map_data);
798

799
	write_record(ncid, NetCDFIO::cf_time, vid_time, 0, number_of_records, dates);
800
	for (map<string, double*>::const_iterator it = map_data.begin(); it != map_data.end(); ++it) {
801
		const string& varname = it->first;
Thomas Egger's avatar
Thomas Egger committed
802
		write_data(ncid, varname, varid[varname], map_data[varname]);
803
804
805
806
807
808
809
810
		delete[] it->second;
	}

	close_file(filename, ncid);

	delete[] dates;
}

811
812
813
// Copy the data from the MeteoData objects into C arrays, perform all necessary
// conversions (multiplications) and set plugin_nodata values where required.
// A loop over all parameters present is performed.
814
815
816
void NetCDFIO::copy_data(const size_t& number_of_stations, const size_t& number_of_records, const std::vector< std::vector<MeteoData> >& vecMeteo,
                         const std::map<size_t, std::string>& map_param_name, std::map<std::string, double*>& map_data_2D)
{
817
	for (map<size_t, string>::const_iterator it = map_param_name.begin(); it != map_param_name.end(); ++it) {
818
819
		const size_t param = it->first;
		const string varname = it->second;
Thomas Egger's avatar
Thomas Egger committed
820
821
822

		bool simple_copy = false, multiply_copy = false;
		double multiplier = IOUtils::nodata;
823
824
825

		double* data = map_data_2D[varname];

Thomas Egger's avatar
Thomas Egger committed
826
827
828
829
830
831
832
833
834
835
		if (param == MeteoData::RH) {
			multiplier = 100.;
			multiply_copy = true;
		} else if (param == MeteoData::HNW) {
			multiply_copy = true;
			multiplier = 1./3600.;
		} else {
			simple_copy = true;
		}

836
837
		for (size_t ii=0; ii<number_of_stations; ++ii) {
			for (size_t jj=0; jj<number_of_records; ++jj) {
Thomas Egger's avatar
Thomas Egger committed
838
				const double& value = vecMeteo[ii][jj](param);
839

Thomas Egger's avatar
Thomas Egger committed
840
				if (value == IOUtils::nodata) {
841
					data[jj*number_of_stations + ii] = plugin_nodata;
Thomas Egger's avatar
Thomas Egger committed
842
843
844
845
846
				} else if (simple_copy) {
					data[jj*number_of_stations + ii] = value;
				} else if (multiply_copy) {
					data[jj*number_of_stations + ii] = value * multiplier;
				}
847
			}
848
849
850
851
		}
	}
}

852
// Create meta data variables in the NetCDF dataset
853
854
void NetCDFIO::create_meta_data(const int& ncid, const int& did, std::map<std::string, double*>& map_data_1D, std::map<std::string, int>& varid)
{
855
	for (map<string, double*>::const_iterator it = map_data_1D.begin(); it != map_data_1D.end(); ++it) {
856
857
858
		int vid;
		const string& varname = it->first;

Thomas Egger's avatar
Thomas Egger committed
859
860
861
862
863
		if (varname == cnrm_timestep) {
			add_0D_variable(ncid, cnrm_timestep, NC_DOUBLE, vid);
		} else {
			add_1D_variable(ncid, varname, NC_DOUBLE, did, vid);
		}
864
		add_attribute(ncid, vid, "_FillValue", plugin_nodata);
Thomas Egger's avatar
Thomas Egger committed
865
		add_attributes_for_variable(ncid, vid, varname);
866
867
868
869
870

		varid[varname] = vid;
	}
}

871
872
// Create the parameter variables in the NetCDF dataset, allocate memory for the
// respective C arrays and store the variable ids in the varid map.
873
874
875
876
void NetCDFIO::create_parameters(const int& ncid, const int& did_time, const int& did_points, const size_t& number_of_records,
						   const size_t& number_of_stations, std::map<size_t, std::string>& map_param_name,
                                 std::map<std::string, double*>& map_data_2D, std::map<std::string, int>& varid)
{
877
	// At this point map_param_name holds all parameters that have values different from nodata
878
	for (map<size_t, string>::iterator it = map_param_name.begin(); it != map_param_name.end();) {
879
		bool create = false;
880
881
		string& varname = it->second;

882
		const map<string, string>::const_iterator it_cnrm = map_name.find(varname);
883
		if (it_cnrm != map_name.end()) {
884
885
886
887
888
889
890
891
892
893
894
895
896
897
			varname = it_cnrm->second; // the offical CNRM name for the parameter
			create = true;
			++it;
		} else {
			if (out_strict) {
				// ignore any parameters not defined in the CNRM standard:
				// if a parameter in map_param_name has no equivalent in the map_name map
				// it is deleted from map_param_name and henceforth ignored.
				map_param_name.erase(it++);
			} else {
				create = true;
				++it;
			}
		}
898

899
		if (create) {
900
901
902
			int vid;

			double* data = new double[number_of_records*number_of_stations];
903
			map_data_2D[varname] = data;
904

905
			add_2D_variable(ncid, varname, NC_DOUBLE, did_time, did_points, vid);
906
			add_attribute(ncid, vid, "_FillValue", plugin_nodata);
907
			add_attributes_for_variable(ncid, vid, varname);
908
909
910

			varid[varname] = vid;
		}
911
	}
912
913
}

914
// Retrieve the parameters in use (parameters, that are different from nodata
915
// for at least one timestamp for at least one station) and store them in
916
// map_param_name. map_param_name associates a MeteoData parameter index with a
917
// string name, that is the CNRM name for the parameter to use in the NetCDF
918
919
920
921
// file. Furthermore this method copies the meta data into the appropriate C
// arrays. The timestep interval is also calculated and added to the map_data_1D
void NetCDFIO::get_parameters(const std::vector< std::vector<MeteoData> >& vecMeteo, std::map<size_t, std::string>& map_param_name,
                              std::map<std::string, double*>& map_data_1D, double*& dates)
922
{
923
	const size_t number_of_records = vecMeteo[0].size();
924
925
	dates = new double[number_of_records];

Thomas Egger's avatar
Thomas Egger committed
926
	double interval = 0;
927
928
	for (size_t ii=0; ii<number_of_records; ii++) {
		dates[ii] = vecMeteo[0][ii].date.getModifiedJulianDate();
929
		if (ii == 1) interval = Optim::round((dates[ii] - dates[ii-1]) * 86400.);
930
931
	}

932
	const size_t nr_of_parameters = (!vecMeteo[0].empty())? vecMeteo[0][0].getNrOfParameters() : 0 ;
933
934
935
936
937
	vector<bool> vec_param_in_use(nr_of_parameters, false);
	vector<string> vec_param_name(nr_of_parameters, "");

	//Check consistency, dates must be existent everywhere
	bool inconsistent = false;
938
	for (size_t ii=0; ii<vecMeteo.size(); ++ii) {
939
		if (number_of_records != vecMeteo[ii].size()) inconsistent = true;
940
		for (size_t jj=0; jj<vecMeteo[ii].size(); ++jj) {
941
942
			const MeteoData& meteo_data = vecMeteo[ii][jj];

943
944
			if (!IOUtils::checkEpsilonEquality(dates[jj], meteo_data.date.getModifiedJulianDate(), NetCDFIO::epsilon))
				inconsistent = true;
945
946

			if (jj == 0) {
947
948
				map_data_1D[cnrm_latitude][ii] = meteo_data.meta.position.getLat();
				map_data_1D[cnrm_longitude][ii] = meteo_data.meta.position.getLon();
Thomas Egger's avatar