WSL/SLF GitLab Repository

NetCDFIO.cc 45.9 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
20
#include <meteoio/ResamplingAlgorithms2D.h>
#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
41
42
43
44
45
46
47
48
 * 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.
 *
 * 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*
49
50
51
52
53
54
55
56
 *
 * @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
57
58
 * - DEMFILE: The filename of the file containing the DEM; [Input] section
 * - DEMVAR: The variable name of the DEM within the DEMFILE; [Input] section
59
 * - METEOFILE: the NetCDF file which shall be used for the meteo parameter input/output; [Input] and [Output] section
60
 * - GRID2DFILE: the NetCDF file which shall be used for gridded input/output; [Input] and [Output] section
61
62
 * - 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
63
64
65
66
67
68
69
70
71
72
73
74
 *
 * @section example Example use
 * @code
 * [Input]
 * DEM     = NETCDF
 * DEMFILE = ./input/Aster_tile.nc
 * DEMVAR  = z
 * @endcode
 *
 * @section Compilation
 * 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.
75
76
 */

77
const double NetCDFIO::plugin_nodata = -9999999.; //CNRM-GAME nodata value
78

79
80
81
82
const std::string NetCDFIO::cf_time = "time";
const std::string NetCDFIO::cf_units = "units";
const std::string NetCDFIO::cf_days = "days since ";
const std::string NetCDFIO::cf_seconds = "seconds since ";
83
84
85
86
87
88
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";
89

90
91
92
const std::string NetCDFIO::cnrm_points = "Number_of_points";
const std::string NetCDFIO::cnrm_latitude = "LAT";
const std::string NetCDFIO::cnrm_longitude = "LON";
93
94
95
const std::string NetCDFIO::cnrm_altitude = "ZS";
const std::string NetCDFIO::cnrm_aspect = "aspect";
const std::string NetCDFIO::cnrm_slope = "slope";
96
97
98
99
100
101
102
103
104
105
106
107
108
109
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
110
const std::string NetCDFIO::cnrm_timestep = "FRC_TIME_STP";
111

112
std::map<std::string, size_t> NetCDFIO::paramname;
113
std::map<std::string, std::string> NetCDFIO::map_name;
114
115
116
117
118
const bool NetCDFIO::__init = NetCDFIO::initStaticData();

bool NetCDFIO::initStaticData()
{
	//Associate unsigned int value and a string representation of a meteo parameter
119
	paramname[cnrm_ta] = MeteoData::TA;
120
121
122
123
	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
124
125
126
127
128
129
130
131
132
	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;
133

134
135
136
137
138
139
	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
140
	map_name["ISWR"] = cnrm_swr_direct;
141
	map_name["HNW"] = cnrm_hnw;
142
143
144
145
	map_name[cnrm_co2air] = cnrm_co2air;
	map_name[cnrm_qair] = cnrm_qair;
	map_name[cnrm_theorsw] = cnrm_theorsw;
	map_name[cnrm_neb] = cnrm_neb;
146

147
148
149
	return true;
}

150
NetCDFIO::NetCDFIO(const std::string& configfile) : cfg(configfile), coordin(), coordinparam(), coordout(), coordoutparam(),
151
                                                    in_dflt_TZ(0.), out_dflt_TZ(0.), in_strict(false), out_strict(false), vecMetaData()
152
153
{
	IOUtils::getProjectionParameters(cfg, coordin, coordinparam, coordout, coordoutparam);
154
	parseInputOutputSection();
155
156
}

157
NetCDFIO::NetCDFIO(const Config& cfgreader) : cfg(cfgreader), coordin(), coordinparam(), coordout(), coordoutparam(),
158
                                              in_dflt_TZ(0.), out_dflt_TZ(0.), in_strict(false), out_strict(false), vecMetaData()
159
160
{
	IOUtils::getProjectionParameters(cfg, coordin, coordinparam, coordout, coordoutparam);
161
	parseInputOutputSection();
162
163
}

164
NetCDFIO::~NetCDFIO() throw() {}
165

166
167
168
169
170
171
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);
172
173
174

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

177
void NetCDFIO::read2DGrid(Grid2DObject& grid_out, const std::string& arguments)
178
{
179
180
181
182
183
184
185
186
	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);
	}
187
188
}

189
void NetCDFIO::read2DGrid(Grid2DObject& grid_out, const MeteoGrids::Parameters& parameter, const Date& date)
190
{
191
	string filename;
192
193
	cfg.getValue("GRID2DFILE", "Input", filename);

194
	const string varname = get_varname(parameter);
195

196
	read2DGrid_internal(grid_out, filename, varname, date);
197
198
}

199
void NetCDFIO::read2DGrid_internal(Grid2DObject& grid_out, const std::string& filename, const std::string& varname, const Date& date)
200
{
201
	const bool is_record = (date != Date());
202
203
	size_t lat_index = 0, lon_index = 1;

204
	int ncid, varid;
205
	vector<int> dimid, dim_varid;
206
207
208
209
210
211
212
	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);

213
	if (is_record) { // In case we're reading a record the first index is always the record index
214
215
216
217
218
219
		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);
	} else if (dimid.size()!=2 || dimlen[lat_index]<2 || dimlen[lon_index]<2) {
220
		throw IOException("Variable '" + varname + "' may only have two dimensions and both have to have length >1", AT);
221
222
223
224
225
	}

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

227
228
	read_data(ncid, dimname[lat_index], dim_varid[lat_index], lat);
	read_data(ncid, dimname[lon_index], dim_varid[lon_index], lon);
229

230
	if (is_record) {
231
		const size_t pos = find_record(ncid, NetCDFIO::cf_time, dimid[0], date.getModifiedJulianDate());
232
233
		if (pos == IOUtils::npos)
			throw IOException("No record for date " + date.toString(Date::ISO), AT);
234

235
236
237
238
		read_data(ncid, varname, varid, pos, dimlen[lat_index], dimlen[lon_index], grid);
	} else {
		read_data(ncid, varname, varid, grid);
	}
239

240
	copy_grid(dimlen[lat_index], dimlen[lon_index], lat, lon, grid, grid_out);
241
242
243

	close_file(filename, ncid);

244
	delete[] lat; delete[] lon; delete[] grid;
245
246
}

247
248
void NetCDFIO::copy_grid(const size_t& latlen, const size_t& lonlen, const double * const lat, const double * const lon,
                         const double * const grid, Grid2DObject& grid_out)
249
{
250
251
	Coords location(coordin, coordinparam);
	location.setLatLon(lat[0], lon[0], grid[0]);
252

253
	double resampling_factor_x = IOUtils::nodata, resampling_factor_y=IOUtils::nodata;
254
	const double cellsize = calculate_cellsize(latlen, lonlen, lat, lon, resampling_factor_x, resampling_factor_y);
255

256
	grid_out.set(lonlen, latlen, cellsize, location);
257

258
259
260
	for (size_t kk=0; kk < latlen; kk++) {
		for (size_t ll=0; ll < lonlen; ll++) {
			grid_out(ll, kk) = IOUtils::standardizeNodata(grid[kk*lonlen + ll], plugin_nodata);
261
262
263
		}
	}

264
265
	if (resampling_factor_x != IOUtils::nodata) {
		grid_out.grid2D = ResamplingAlgorithms2D::BilinearResampling(grid_out.grid2D, resampling_factor_x, resampling_factor_y);
266
267
268
		grid_out.ncols = grid_out.grid2D.getNx();
		grid_out.nrows = grid_out.grid2D.getNy();
	}
269
270
}

271
272
273
/* 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
274
275
 * 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
276
277
 * quadratic cells.
 */
278
279
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)
280
{
281
282
	Coords llcorner(coordin, coordinparam);
	llcorner.setLatLon(lat[0], lon[0], IOUtils::nodata);
283
284
285
286

	Coords urcorner(coordin, coordinparam);
	urcorner.setLatLon(lat[latlen-1], lon[lonlen-1], IOUtils::nodata);

287
288
289
290
291
292
293
294
295
296
297
298
	const double ll_easting=llcorner.getEasting(), ll_northing=llcorner.getNorthing();
	const double ur_easting=urcorner.getEasting(), ur_northing=urcorner.getNorthing();

	const double distanceX = ur_easting - ll_easting;
	const double distanceY = ur_northing - ll_northing;
	if(distanceX<0 || distanceY<0) {
		ostringstream ss;
		ss << "Can not compute cellsize: this is most probably due to an inappropriate input coordinate system (COORDSYS).";
		ss << "Please configure one that can accomodate (" << llcorner.getLat() << "," << llcorner.getLon() << ") - ";
		ss << "(" << urcorner.getLat() << "," << urcorner.getLon() << ")";
		throw InvalidArgumentException(ss.str(), AT);
	}
299

300
	// lonlen, latlen are decremented by 1; n linearly connected points have (n-1) connections
301
302
	const double cellsize_x = distanceX / (lonlen-1);
	const double cellsize_y = distanceY / (latlen-1);
303

304
	// We're using a precision for the cellsize that is equal to 1cm, more
305
	// precision makes ensuing calculations numerically instable
306
307
	const double value =  max(cellsize_x, cellsize_y) * 100.0;
	const double cellsize = floor(value) / 100.0;
308

309
310
	if (cellsize_x == cellsize_y) {
		return cellsize_x;
311
	} else {
312
313
314
315
		factor_x =  cellsize_x / cellsize;
		factor_y =  cellsize_y / cellsize;

		return cellsize;
316
	}
317
318
}

319
320
void NetCDFIO::readDEM(DEMObject& dem_out)
{
321
	string filename, varname;
322

323
324
	cfg.getValue("DEMFILE", "Input", filename);
	cfg.getValue("DEMVAR", "Input", varname);
325

326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
	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)
{
343
	if (!vecMetaData.empty()) { // We already have meta data
344
345
346
347
		vecStation = vecMetaData;
		return;
	}

348
	string filename;
349
350
351
352
353
354
355
356
357
358
359
360
361
	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)
{
362
363
364
	vecStation.clear();

	int dimid;
365
	size_t dimlen;
366
	map<string, int> map_vid;
367

368
	get_dimension(ncid, cnrm_points, dimid, dimlen);
369
370
371
	if (dimlen == 0) return; // There are no stations

	get_meta_data_ids(ncid, map_vid);
372
373
374
375

	double *alt = new double[dimlen];
	double *lat = new double[dimlen];
	double *lon = new double[dimlen];
376
377
	double *aspect = new double[dimlen];
	double *slope = new double[dimlen];
378

379
380
381
382
383
	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);
384
385
386
387
388
389

	//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]);
390

391
		ss << (ii+1);
392
		const string id( ss.str() );
393
394
395
		ss.str("");

		ss << "Station " << (ii +1);
396
		const string name( ss.str() );
397
		ss.str("");
398

399
		StationData tmp(location, id, name);
400
		const double aspect_bearing = (aspect[ii] < 0) ? 0 : aspect[ii]; // aspect allowed to be -1 in CNRM format...
401
		tmp.setSlope(slope[ii], aspect_bearing);
402
403
404
		vecStation.push_back(tmp);
	}

405
	delete[] alt; delete[] lat; delete[] lon; delete[] aspect; delete[] slope;
406
407
}

408
409
410
411
412
413
414
415
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

416
	for (vector<string>::const_iterator it = varname.begin(); it != varname.end(); ++it) {
417
418
419
420
		int varid;
		const string& name = *it;

		get_variable(ncid, name, varid);
421
422
		if (!check_dimensions(ncid, name, varid, dimensions))
			throw IOException("Variable '" + name  + "' fails dimension check", AT);
423
424
425
426
427

		map_vid[name] = varid;
	}
}

428
void NetCDFIO::readMeteoData(const Date& dateStart, const Date& dateEnd, std::vector< std::vector<MeteoData> >& vecMeteo, const size_t&)
429
{
430
431
	vecMeteo.clear();

432
	string filename;
433
434
435
436
437
438
439
440
441
442
443
	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
444
445

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

450
			readData(ncid, index_start, vec_date, map_parameters, meteo_data, vecMeteo);
451
452
453
454
455
456
		}
	}

	close_file(filename, ncid);
}

457
void NetCDFIO::readData(const int& ncid, const size_t& index_start, const std::vector<Date>& vec_date,
458
                        const std::map<std::string, size_t>& map_parameters, const MeteoData& meteo_data, std::vector< std::vector<MeteoData> >& vecMeteo)
459
{
460
461
	const size_t number_of_stations = vecMetaData.size();
	const size_t number_of_records = vec_date.size();
462

463
	// Allocate all the MeteoData objects based on the template meteo_data
464
	vector<MeteoData> tmp_vec(number_of_records, meteo_data);
465
466
467
468
469
470
471
	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);
	}

472
	// Allocate enough linear space for each parameter and read the data from NetCDF
473
	map<string, double*> map_data;
474
	for (map<string, size_t>::const_iterator it = map_parameters.begin(); it != map_parameters.end(); ++it) {
475
		double* data = new double[number_of_stations*number_of_records];
476
		const string& varname = it->first;
477

478
		map_data[varname] = data;
479

480
		int varid;
481
482
483
484
		get_variable(ncid, varname, varid);
		read_data_2D(ncid, varname, varid, index_start, number_of_records, number_of_stations, data);
	}

485
	copy_data(ncid, map_parameters, map_data, number_of_stations, number_of_records, vecMeteo);
486

487
	for (map<string, double*>::const_iterator it = map_data.begin(); it != map_data.end(); ++it) {
488
489
490
491
		delete[] it->second;
	}
}

492
493
494
495
496
497
// 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.
498
void NetCDFIO::copy_data(const int& ncid, const std::map<std::string, size_t>& map_parameters, const std::map<std::string, double*> map_data,
499
                         const size_t& number_of_stations, const size_t& number_of_records, std::vector< std::vector<MeteoData> >& vecMeteo)
500
{
501
	for (map<string, double*>::const_iterator it = map_data.begin(); it != map_data.end(); ++it) {
502
503
504
		const string& varname = it->first;

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

		if (param == IOUtils::npos) {
510
511
			if ((varname == cnrm_snowf) || (varname == cnrm_hnw)) {
				int varid;
Thomas Egger's avatar
Thomas Egger committed
512
513
				get_variable(ncid, cnrm_timestep, varid);
				read_value(ncid, cnrm_timestep, varid, multiplier);
514

Thomas Egger's avatar
Thomas Egger committed
515
				if (multiplier <= 0) throw InvalidArgumentException("The variable '" + cnrm_timestep + "' is invalid", AT);
516
517
518
519

				hnw_measurement = true;
			} else if ((varname == cnrm_swr_diffuse) || (varname == cnrm_swr_direct)) {
				sw_measurement = true;
520
521
			} else {
				throw IOException("Don't know how to deal with parameter " + varname, AT);
522
523
			}
		} else {
524
			if (varname == cnrm_rh) {
525
526
527
528
529
530
531
				mutiply_copy = true;
				multiplier = 0.01;
			} else {
				simple_copy = true;
			}
		}

532
		// Loop through all times and all stations
533
		for (size_t jj=0; jj<number_of_records; jj++) {
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
			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;
					}
551
				} else if (hnw_measurement) {
552
					if (!nodata) {
553
						double& hnw = vecMeteo[ii][jj](MeteoData::HNW);
554
555
556
						if (hnw == IOUtils::nodata) hnw = 0.0;
						hnw += value * multiplier;
					}
557
558
559
560
				} else if (sw_measurement) {
					if (!nodata) {
						double& iswr = vecMeteo[ii][jj](MeteoData::ISWR);
						if (iswr == IOUtils::nodata) iswr = 0.0;
561
						iswr += value;
562
					}
563
				}
564
565
566
567
568
			}
		}
	}
}

569
570
571
// 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
572
573
// 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().
574
void NetCDFIO::get_parameters(const int& ncid, std::map<std::string, size_t>& map_parameters, MeteoData& meteo_data)
575
{
576
577
	vector<string> dimensions;
	dimensions.push_back(cf_time);
578
579
580
581
582
583
584
	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;
585
		//cout << "Found parameter: " << name << endl;
586

587
588
589
590
		// Check if parameter exists in paramname, which holds strict CNRM parameters
		map<string, size_t>::const_iterator strict_it = paramname.find(name);
		if (strict_it != paramname.end()) { // parameter is a part of the CNRM specification
			size_t index = strict_it->second;
591
592
593
594
595

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

596
597
598
			map_parameters[name] = index;
		} else if (!in_strict) { // parameter will be read anyway
			size_t index = IOUtils::npos;
599

600
601
602
603
604
605
606
			if (meteo_data.param_exists(name)) {
				index = meteo_data.getParameterIndex(name);
			} else {
			 	index = meteo_data.addParameter(name);
			}

			map_parameters[name] = index;
607
608
609
610
		}
	}
}

611
612
613
614
// 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
615
616
void NetCDFIO::get_indices(const int& ncid, const Date& dateStart, const Date& dateEnd, size_t& indexStart, size_t& indexEnd, std::vector<Date>& vecDate)
{
617
618
	indexStart = indexEnd = IOUtils::npos;

619
620
621
622
623
	int varid, dimid;
	size_t dimlen;
	get_dimension(ncid, NetCDFIO::cf_time, dimid, dimlen);
	get_variable(ncid, NetCDFIO::cf_time, varid);

624
	// Get the units attribute and calculate the offset date
625
626
627
	string units_str;
	NetCDFIO::TimeUnit unit_type;
	Date offset;
628
629
	get_attribute(ncid, NetCDFIO::cf_time, varid, cf_units, units_str);
	calculate_offset(units_str, unit_type, offset);
630
631
632
633

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

634
	// Firstly, check whether search makes any sense, that is dateStart and dateEnd overlap with the times present
635
636
637
	bool search = true;
	if (dimlen > 0) {
		Date time_start(offset), time_end(offset);
638

639
640
641
642
643
644
645
		double start = time[0];
		double end = time[dimlen-1];

		if (unit_type == seconds) {
			start /= 86400;
			end   /= 86400;
		}
646
647
		time_start += Date(start, 0.0);
		time_end += Date(end, 0.0);
648
649
650
651
652

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

653
654
	// If search is feasible then loop through the existent timestamps and find the relevant indices
	bool start_found = false;
655
656
657
658
659
660
	if (search) {
		for (size_t ii=0; ii<dimlen; ii++) {
			if (unit_type == seconds) {
				time[ii] /= 86400;
			}

661
			const Date tmp_date = offset + Date(time[ii], 0.0);
662

663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
			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;
		}
	}
678

679
680
681
	delete[] time;
}

682
683
// The CNRM timestamps have an offset that is saved in the units attribute of
// the time variable - this method retrieves that offset
684
685
686
void NetCDFIO::calculate_offset(const std::string& units, NetCDFIO::TimeUnit& time_unit, Date& offset)
{
	string tmp(units);
687
688
	const size_t found_sec = units.find(NetCDFIO::cf_seconds);
	const size_t found_day = units.find(NetCDFIO::cf_days);
689
690
691
692
693
694
695
696

	if (found_sec != string::npos) {
		time_unit = seconds;
		tmp = tmp.substr(found_sec + NetCDFIO::cf_seconds.size());
	} else if (found_day != string::npos) {
		time_unit = days;
		tmp = tmp.substr(found_day+ + NetCDFIO::cf_days.size());
	} else {
697
		throw InvalidFormatException("Variable '"+NetCDFIO::cf_time+"' has no valid attribute '" + cf_units + "'" , AT);
698
	}
699

700
	const bool success = IOUtils::convertString(offset, tmp, in_dflt_TZ);
701
	if (!success) throw InvalidFormatException("Cannot parse time: " + tmp, AT);
702
703
}

704
void NetCDFIO::writeMeteoData(const std::vector< std::vector<MeteoData> >& vecMeteo, const std::string&)
705
{
706
	const size_t number_of_stations = vecMeteo.size();
707
708
	if (number_of_stations == 0) return; //Nothing to write

709
	const size_t number_of_records = vecMeteo[0].size();
710

711
	string filename;
712
713
714
715
	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;
716

717
	const bool exists = IOUtils::fileExists(filename);
718
	if (exists) remove(filename.c_str()); // NOTE: file is deleted if it exists
719
720

	double* dates;
721
722
723
	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
724
725
726
	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];
727
728
729
730

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

Thomas Egger's avatar
Thomas Egger committed
731
	get_parameters(vecMeteo, map_param_name, map_data, dates);
732

733
734
	create_file(filename, NC_CLASSIC_MODEL, ncid);
	create_time = create_points = create_locations = create_variables = true;
735
736
737

	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
738
739
	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);
740
741
742

	end_definitions(filename, ncid);

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

745
	write_record(ncid, NetCDFIO::cf_time, vid_time, 0, number_of_records, dates);
746
	for (map<string, double*>::const_iterator it = map_data.begin(); it != map_data.end(); ++it) {
747
		const string& varname = it->first;
Thomas Egger's avatar
Thomas Egger committed
748
		write_data(ncid, varname, varid[varname], map_data[varname]);
749
750
751
752
753
754
755
756
		delete[] it->second;
	}

	close_file(filename, ncid);

	delete[] dates;
}

757
758
759
// 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.
760
761
762
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)
{
763
	for (map<size_t, string>::const_iterator it = map_param_name.begin(); it != map_param_name.end(); ++it) {
764
765
		const size_t& param = it->first;
		const string& varname = it->second;
Thomas Egger's avatar
Thomas Egger committed
766
767
768

		bool simple_copy = false, multiply_copy = false;
		double multiplier = IOUtils::nodata;
769
770
771

		double* data = map_data_2D[varname];

Thomas Egger's avatar
Thomas Egger committed
772
773
774
775
776
777
778
779
780
781
		if (param == MeteoData::RH) {
			multiplier = 100.;
			multiply_copy = true;
		} else if (param == MeteoData::HNW) {
			multiply_copy = true;
			multiplier = 1./3600.;
		} else {
			simple_copy = true;
		}

782
783
		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
784
				const double& value = vecMeteo[ii][jj](param);
785

Thomas Egger's avatar
Thomas Egger committed
786
				if (value == IOUtils::nodata) {
787
					data[jj*number_of_stations + ii] = plugin_nodata;
Thomas Egger's avatar
Thomas Egger committed
788
789
790
791
792
				} else if (simple_copy) {
					data[jj*number_of_stations + ii] = value;
				} else if (multiply_copy) {
					data[jj*number_of_stations + ii] = value * multiplier;
				}
793
			}
794
795
796
797
		}
	}
}

798
// Create meta data variables in the NetCDF dataset
799
800
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)
{
801
	for (map<string, double*>::const_iterator it = map_data_1D.begin(); it != map_data_1D.end(); ++it) {
802
803
804
		int vid;
		const string& varname = it->first;

Thomas Egger's avatar
Thomas Egger committed
805
806
807
808
809
		if (varname == cnrm_timestep) {
			add_0D_variable(ncid, cnrm_timestep, NC_DOUBLE, vid);
		} else {
			add_1D_variable(ncid, varname, NC_DOUBLE, did, vid);
		}
810
		add_attribute(ncid, vid, "_FillValue", plugin_nodata);
Thomas Egger's avatar
Thomas Egger committed
811
		add_attributes_for_variable(ncid, vid, varname);
812
813
814
815
816

		varid[varname] = vid;
	}
}

817
818
// Create the parameter variables in the NetCDF dataset, allocate memory for the
// respective C arrays and store the variable ids in the varid map.
819
820
821
822
823
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)
{
	map<string, string>::const_iterator it_cnrm;
824

825
	// At this point map_param_name holds all parameters that have values different from nodata
826
	for (map<size_t, string>::iterator it = map_param_name.begin(); it != map_param_name.end();) {
827
		bool create = false;
828
829
830
831
		string& varname = it->second;

		it_cnrm = map_name.find(varname);
		if (it_cnrm != map_name.end()) {
832
833
834
835
836
837
838
839
840
841
842
843
844
845
			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;
			}
		}
846

847
		if (create) {
848
849
850
			int vid;

			double* data = new double[number_of_records*number_of_stations];
851
			map_data_2D[varname] = data;
852

853
			add_2D_variable(ncid, varname, NC_DOUBLE, did_time, did_points, vid);
854
			add_attribute(ncid, vid, "_FillValue", plugin_nodata);
855
			add_attributes_for_variable(ncid, vid, varname);
856
857
858

			varid[varname] = vid;
		}
859
	}
860
861
}

862
// Retrieve the parameters in use (parameters, that are different from nodata
863
// for at least one timestamp for at least one station) and store them in
864
// map_param_name. map_param_name associates a MeteoData parameter index with a
865
// string name, that is the CNRM name for the parameter to use in the NetCDF
866
867
868
869
// 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)
870
{
871
	const size_t number_of_records = vecMeteo[0].size();
872
873
	dates = new double[number_of_records];

Thomas Egger's avatar
Thomas Egger committed
874
	double interval = 0;
875
876
	for (size_t ii=0; ii<number_of_records; ii++) {
		dates[ii] = vecMeteo[0][ii].date.getModifiedJulianDate();
877
		if (ii == 1) interval = Optim::round((dates[ii] - dates[ii-1]) * 86400.);
878
879
880
	}

	size_t nr_of_parameters = 0;
881
	if (!vecMeteo[0].empty()) nr_of_parameters = vecMeteo[0][0].getNrOfParameters();
882
883
884
885
886
887

	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;
888
	for (size_t ii=0; ii<vecMeteo.size(); ++ii) {
889
		if (number_of_records != vecMeteo[ii].size()) inconsistent = true;
890
		for (size_t jj=0; jj<vecMeteo[ii].size(); ++jj) {
891
892
893
894
895
			const MeteoData& meteo_data = vecMeteo[ii][jj];

			if (dates[jj] != meteo_data.date.getModifiedJulianDate()) inconsistent = true;

			if (jj == 0) {
896
897
				map_data_1D[cnrm_latitude][ii] = meteo_data.meta.position.getLat();
				map_data_1D[cnrm_longitude][ii] = meteo_data.meta.position.getLon();
898
899
900
901
902
903
				map_data_1D[cnrm_altitude][ii] = meteo_data.meta.position.getAltitude();
				map_data_1D[cnrm_slope][ii] = meteo_data.meta.getSlopeAngle();
				map_data_1D[cnrm_aspect][ii] = meteo_data.meta.getAzimuth();
			}

			//Check which parameters are in use
904
			for (size_t kk=0; kk<nr_of_parameters; ++kk) {
905
906
907
908
909
910
911
912
913
914
915
916
				if (!vec_param_in_use[kk]){
					if (meteo_data(kk) != IOUtils::nodata){
						vec_param_in_use[kk] = true;
						vec_param_name[kk] = meteo_data.getNameForParameter(kk);
					}
				}
			}
		}
	}

	if (inconsistent) throw IOException("Inconsistent dates in vecMeteo between different stations", AT);

917
	for (size_t kk=0; kk<nr_of_parameters; ++kk) {
918
919
920
		if (vec_param_in_use[kk])
			map_param_name[kk] = vec_param_name[kk];
	}
Thomas Egger's avatar
Thomas Egger committed
921
922
923
924
925


	double* timestep = new double[1];
	*timestep = interval;
	map_data_1D[cnrm_timestep] = timestep;
926
927
928
929
930
931
932
933
934
935
}

void NetCDFIO::readPOI(std::vector<Coords>&)
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

void NetCDFIO::write2DGrid(const Grid2DObject& grid_in, const std::string& arguments)
{
936
	// arguments is a string of the format filname:varname
937
938
939
940
941
942
	vector<string> vec_argument;
	IOUtils::readLineToVec(arguments, vec_argument, ':');

	if (vec_argument.size() != 2)
		throw InvalidArgumentException("The format for the arguments to NetCDFIO::write2DGrid is filename:varname", AT);

943
	write2DGrid_internal(grid_in, vec_argument[0], vec_argument[1]);
944
945
946
947
}

void NetCDFIO::write2DGrid(const Grid2DObject& grid_in, const MeteoGrids::Parameters& parameter, const Date& date)
{
948
	string filename;
949
950
	cfg.getValue("GRID2DFILE", "Output", filename);

951
	const string varname = get_varname(parameter);
952
953
954
955
956
957

	write2DGrid_internal(grid_in, filename, varname, date);
}

void NetCDFIO::write2DGrid_internal(const Grid2DObject& grid_in, const std::string& filename, const std::string& varname, const Date& date)
{
958
959
	const bool is_record = (date != Date());
	const bool exists = IOUtils::fileExists(filename);
960

961
962
963
964
965
966
967
968
969
970
971
972
	double *lat_array = new double[grid_in.nrows];
	double *lon_array = new double[grid_in.ncols];
	double *data = new double[grid_in.nrows * grid_in.ncols];

	calculate_dimensions(grid_in, lat_array, lon_array);
	fill_data(grid_in, data);

	int ncid, did_lat, did_lon, did_time, vid_lat, vid_lon, vid_var, vid_time;
	bool create_dimensions(false), create_variable(false), create_time(false);

	if (exists) {
		open_file(filename, NC_WRITE, ncid);
973

974
		//check of lat/lon are defined and consistent