WSL/SLF GitLab Repository

GRIBIO.cc 41.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/***********************************************************************************/
/*  Copyright 2012 WSL Institute for Snow and Avalanche Research    SLF-DAVOS      */
/***********************************************************************************/
/* 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 "GRIBIO.h"

20
21
#include <meteoio/meteolaws/Atmosphere.h>
#include <meteoio/meteolaws/Meteoconst.h> //for PI
22
#include <meteoio/DEMObject.h>
23

24
#include <cmath>
25
26
27
28
29
30
31
#include <errno.h>
#include <grib_api.h>

using namespace std;

namespace mio {
/**
32
 * @page gribio GRIBIO
33
 * @section gribio_format Format and limitations
34
 * This plugin reads GRIB (https://en.wikipedia.org/wiki/GRIB) files as produced by meteorological models.
35
 * Being based on GRIB API (http://www.ecmwf.int/products/data/software/grib_api.html), it should support both version 1 and 2 of the format (please note that grib_api must be compiled with Position Independent Code ("fPIC" flag)).
36
 * Fields are read based on their marsParam code (this is built as {grib parameter number}.{grib table number} the table being preferably table 2, the parameter being preferably WMO standardized, as in http://dss.ucar.edu/docs/formats/grib/gribdoc/params.html) and levels
37
 * (levels description is available at http://www.nco.ncep.noaa.gov/pmb/docs/on388/).
38
 *
39
40
41
42
43
44
45
46
47
 * Several assumptions/approximations are held/made when reading grids:
 * - since models usually use rotated latitude/longitude (see http://www.cosmo-model.org/content/model/documentation/core/default.htm, part I, chapter 3.3), the center of the domain can be approximated by a tangential cartesian coordinate system. We therefore don't re-project the lat/lon grid and use it "as is".
 * - however, no correction for grid rotation is currently performed. If a grid rotation is specified on top of the rotated coordinate system, an error message will be given
 * - the cell size is computed at the center of the domain. This is performed by retrieving the latitude and longitude increments in the rotated coordinates, computing the point at center+increment in geographic coordinates and computing the equivalent geographic latitude and longitude increment. These increments are then converted to distances along the parallel and meridian at the true center latitude (see https://en.wikipedia.org/wiki/Latitude#The_length_of_a_degree_of_latitude).
 * - the average cell size (ie: average between x and y) is used to move the center point to the lower left corner. This will be returned as lower left corner geolocalization of the grid.
 *
 * This means that close to the center of the grid, coordinates and distances will work as expected, but the distortion will increase when moving away from the center and can become significant. As examples for domain size, cone can look at the MeteoSwiss domain definition at http://www.cosmo-model.org/content/tasks/operational/meteoSwiss/default.htm.
 *
 * As a side note, when calling read2DGrid(grid, filename), it will returns the first grid that is found.
48
 *
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
 * @section cosmo_partners COSMO Group
 * This plugin has been developed primarily for reading GRIB files produced by COSMO (http://www.cosmo-model.org/) at MeteoSwiss.
 * COSMO (COnsortium for Small scale MOdelling) represents a non-hydrostatic limited-area atmospheric model, to be used both for operational and for research applications by the members of the consortium. The Consortium has the following members:
 *  - Germany, DWD, Deutscher Wetterdienst
 *  - Switzerland, MCH, MeteoSchweiz
 *  - Italy, USAM, Ufficio Generale Spazio Aereo e Meteorologia
 *  - Greece, HNMS, Hellenic National Meteorological Service
 *  - Poland, IMGW, Institute of Meteorology and Water Management
 *  - Romania, NMA, National Meteorological Administration
 *  - Russia, RHM, Federal Service for Hydrometeorology and Environmental Monitoring
 *  - Germany, AGeoBw, Amt für GeoInformationswesen der Bundeswehr
 *  - Italy, CIRA, Centro Italiano Ricerche Aerospaziali
 *  - Italy, ARPA-SIMC, ARPA Emilia Romagna Servizio Idro Meteo Clima
 *  - Italy, ARPA Piemonte, Agenzia Regionale per la Protezione Ambientale Piemonte
 *
64
65
 * @section gribio_units Units
 * As specified by WMO.
66
 *
67
 * @section gribio_keywords Keywords
68
 * This plugin uses the following keywords:
69
70
71
72
73
 * - COORDSYS: coordinate system (see Coords)
 * - COORDPARAM: extra coordinates parameters (see Coords)
 * - GRID2DPATH: path where to find the grids
 * - GRIB_DEM_UPDATE: recompute slope/azimuth from the elevations when reading a DEM (default=false,
 * that is we use the slope and azimuth included in the GRIB file)
74
 * - GRIB_PREFIX: prefix to append when generating a file name for reading (ie: something like "laf" for Cosmo-Analysis-full domain), optional
75
 * - GRIB_EXT: grib file extension, or <i>none</i> for no file extension (default: .grb)
76
 * - STATION#: coordinates for virtual stations (if using GRIB as METEO plugin). Each station is given by its coordinates and the closest
77
78
 * grid point will be chosen. Coordinates are given one one line as "lat lon" or "xcoord ycoord epsg_code". If a point leads to duplicate grid points,
 * it will be removed from the list.
79
 *
80
81
82
 */

const double GRIBIO::plugin_nodata = -999.; //plugin specific nodata value. It can also be read by the plugin (depending on what is appropriate)
83
const double GRIBIO::tz_in = 0.; //GRIB time zone, always UTC
84
const std::string GRIBIO::default_ext=".grb"; //filename extension
85
86
const double GRIBIO::to_rad = Cst::PI / 180.0;
const double GRIBIO::to_deg = 180.0 / Cst::PI;
87
88
89
90

GRIBIO::GRIBIO(void (*delObj)(void*), const Config& i_cfg) : IOInterface(delObj), cfg(i_cfg)
{
	setOptions();
91
92
93
	indexed = false;
	idx=NULL;
	fp = NULL;
94
	meteo_initialized = false;
95
	latitudeOfNorthernPole = longitudeOfNorthernPole = IOUtils::nodata;
96
97
	bearing_offset = IOUtils::nodata;
	cellsize_x = cellsize_y = IOUtils::nodata;
98
99
100
101
102
}

GRIBIO::GRIBIO(const std::string& configfile) : IOInterface(NULL), cfg(configfile)
{
	setOptions();
103
104
105
	indexed = false;
	idx=NULL;
	fp = NULL;
106
	meteo_initialized = false;
107
	latitudeOfNorthernPole = longitudeOfNorthernPole = IOUtils::nodata;
108
109
	bearing_offset = IOUtils::nodata;
	cellsize_x = cellsize_y = IOUtils::nodata;
110
111
112
113
114
}

GRIBIO::GRIBIO(const Config& cfgreader) : IOInterface(NULL), cfg(cfgreader)
{
	setOptions();
115
116
117
	indexed = false;
	idx=NULL;
	fp = NULL;
118
	meteo_initialized = false;
119
	latitudeOfNorthernPole = longitudeOfNorthernPole = IOUtils::nodata;
120
121
	bearing_offset = IOUtils::nodata;
	cellsize_x = cellsize_y = IOUtils::nodata;
122
123
124
125
}

GRIBIO::~GRIBIO() throw()
{
126
	cleanup();
127
128
129
130
}

void GRIBIO::setOptions()
{
131
	std::string coordout, coordoutparam;
132
	IOUtils::getProjectionParameters(cfg, coordin, coordinparam, coordout, coordoutparam);
133
	update_dem = false;
134
135
136
137
138
139
140

	string tmp="";
	cfg.getValue("GRID2D", "Input", tmp, Config::nothrow);
	if (tmp == "GRIB") { //keep it synchronized with IOHandler.cc for plugin mapping!!
		cfg.getValue("GRID2DPATH", "Input", grid2dpath_in);
		cfg.getValue("GRIB_DEM_UPDATE", "Input", update_dem, Config::nothrow);
	}
141
	cfg.getValue("GRIB_PREFIX", "Input", prefix, Config::nothrow);
142
143

	ext = default_ext;
144
	cfg.getValue("GRIB_EXT", "Input", ext, Config::nothrow);
145
	if(ext=="none") ext="";
146
}
147

148
149
150
151
152
153
154
155
156
157
void GRIBIO::readStations()
{
	cfg.getValue("METEOPATH", "Input", meteopath_in);
	size_t current_stationnr = 1;
	string current_station;
	do {
		current_station = "";
		stringstream ss;
		ss << "STATION" << current_stationnr;
		cfg.getValue(ss.str(), "Input", current_station, Config::nothrow);
158
		IOUtils::stripComments(current_station);
159
160
161
162
163
164
165

		if (current_station != "") {
			addStation(current_station);
			std::cout <<  "\tRead virtual station " << vecPts.back().printLatLon() << "\n";
		}
		current_stationnr++;
	} while (current_station != "");
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
}

void GRIBIO::addStation(const std::string& coord_spec)
{
	std::istringstream iss(coord_spec);
	double coord1=IOUtils::nodata, coord2=IOUtils::nodata;
	int epsg=IOUtils::inodata;

	iss >> std::skipws >> coord1;
	iss >> std::skipws >> coord2;
	iss >> std::skipws >> epsg;

	if(coord1!=IOUtils::nodata && coord2!=IOUtils::nodata && epsg!=IOUtils::inodata) {
		Coords point;
		point.setEPSG(epsg);
		point.setXY(coord1, coord2, IOUtils::nodata);
		vecPts.push_back(point);
		return;
	}
	if(coord1!=IOUtils::nodata && coord2!=IOUtils::nodata) {
		Coords point(coordin, coordinparam);
		point.setLatLon(coord1, coord2, IOUtils::nodata);
		vecPts.push_back(point);
		return;
	}

	throw InvalidArgumentException("Coordinate specification \""+coord_spec+"\" is invalid!", AT);
193
194
}

195
void GRIBIO::listKeys(grib_handle** h, const std::string& filename)
196
{
197
	const unsigned int MAX_VAL_LEN=1024; //max value string length in GRIB
198
	unsigned long key_iterator_filter_flags=GRIB_KEYS_ITERATOR_ALL_KEYS;
199
200
201
202
203
204
205
	char* name_space=NULL;
	grib_keys_iterator* kiter=grib_keys_iterator_new(*h,key_iterator_filter_flags,name_space);

	if (!kiter) {
		cleanup();
		throw IOException("Unable to create keys iterator for \""+filename+"\"", AT);
	}
206

207
208
209
210
211
212
213
214
215
216
217
218
219
220
	//Iterating over all keys
	while(grib_keys_iterator_next(kiter)) {
		char value[MAX_VAL_LEN];
		size_t vlen=MAX_VAL_LEN;
		const char* name = grib_keys_iterator_get_name(kiter);
		vlen=MAX_VAL_LEN;
		bzero(value,vlen);
		GRIB_CHECK(grib_get_string(*h,name,value,&vlen),name);
		std::cout << name << " = " << value << "\n";
	}

	grib_keys_iterator_delete(kiter);
}

221
void  GRIBIO::listFields(const std::string& filename)
222
{
223
	grib_handle* h=NULL;
224
	int err=0;
225

226
	//loop over the messages
227
228
229
230
231
232
	while((h = grib_handle_new_from_file(0,fp,&err)) != NULL) {
		if(!h) {
			cleanup();
			throw IOException("Unable to create grib handle for \""+filename+"\"", AT);
		}

233
234
235
236
237
238
239
240
		double marsParam;
		GRIB_CHECK(grib_get_double(h,"marsParam",&marsParam),0);

		size_t len=500;
		char name[500], shortname[500];
		GRIB_CHECK(grib_get_string(h,"name",name, &len),0);
		len=500; //this has been reset by the previous call... no comments
		GRIB_CHECK(grib_get_string(h,"shortName",shortname, &len),0);
241
		long levelType;
242
243
244
245
		GRIB_CHECK(grib_get_long(h,"indicatorOfTypeOfLevel", &levelType),0);
		long level=0;
		if(levelType!=1) GRIB_CHECK(grib_get_long(h,"level", &level),0);
		std::cout << marsParam << " " << shortname << " " << name << " type " << levelType << " level " << level << "\n";
246
		grib_handle_delete(h);
247
248
	}
}
249

250
void GRIBIO::getDate(grib_handle* h, Date &base, double &d1, double &d2) {
251
252
253
254
255
256
	long dataDate, dataTime;
	GRIB_CHECK(grib_get_long(h,"dataDate",&dataDate),0);
	GRIB_CHECK(grib_get_long(h,"dataTime",&dataTime),0);

	const int year=dataDate/10000, month=dataDate/100-year*100, day=dataDate-month*100-year*10000;
	const int hour=dataTime/100, minutes=dataTime-hour*100;
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
	base.setDate(year, month, day, hour, minutes, tz_in);

	//reading offset to base date/time, as used for forecast, computed at time t for t+offset
	long stepUnits, startStep, endStep;
	GRIB_CHECK(grib_get_long(h,"stepUnits",&stepUnits),0);
	GRIB_CHECK(grib_get_long(h,"startStep",&startStep),0);
	GRIB_CHECK(grib_get_long(h,"endStep",&endStep),0);

	double step_units; //in julian, ie. in days
	switch(stepUnits) {
		case 0: //minutes
			step_units = 1./(24.*60.);
			break;
		case 1: //hours
			step_units = 1./24.;
			break;
		case 2: //days
			step_units = 1.;
			break;
		case 10: //3 hours
			step_units = 3./24.;
			break;
		case 11: //6 hours
			step_units = 6./24.;
			break;
		case 12: //12 hours
			step_units = 12./24.;
			break;
		case 13: //seconds
			step_units = 1./(24.*3600.);
			break;
		case 14: //15 minutes
			step_units = 15./(24.*60.);
			break;
		case 15: //30 minutes
			step_units = 30./(24.*60.);
			break;
		default:
			std::stringstream ss;
			ss << "GRIB file using stepUnits=" << stepUnits << ", which is not supported";
			throw InvalidFormatException(ss.str(), AT);
	}

	d1 = startStep*step_units;
	d2 = endStep*step_units;
302
303
}

304
Coords GRIBIO::getGeolocalization(grib_handle* h, double &cell_x, double &cell_y)
305
{
306
	//getting transformation parameters
307
308
309
310
	double angleOfRotationInDegrees;
	GRIB_CHECK(grib_get_double(h,"angleOfRotationInDegrees",&angleOfRotationInDegrees),0);
	if(angleOfRotationInDegrees!=0.) {
		throw InvalidArgumentException("Rotated grids not supported!", AT);
311
	}
312
	double latitudeOfSouthernPole, longitudeOfSouthernPole;
313
314
	GRIB_CHECK(grib_get_double(h,"latitudeOfSouthernPoleInDegrees",&latitudeOfSouthernPole),0);
	GRIB_CHECK(grib_get_double(h,"longitudeOfSouthernPoleInDegrees",&longitudeOfSouthernPole),0);
315
316
	latitudeOfNorthernPole = -latitudeOfSouthernPole;
	longitudeOfNorthernPole = longitudeOfSouthernPole+180.;
317

318
319
	//determining llcorner, urcorner and center coordinates
	double ll_latitude, ll_longitude, ll_lat, ll_lon;
320
321
	GRIB_CHECK(grib_get_double(h,"latitudeOfFirstGridPointInDegrees",&ll_latitude),0);
	GRIB_CHECK(grib_get_double(h,"longitudeOfFirstGridPointInDegrees",&ll_longitude),0);
322
323
	Coords::rotatedToTrueLatLon(latitudeOfNorthernPole, longitudeOfNorthernPole, ll_latitude, ll_longitude, ll_lat, ll_lon);
	double ur_latitude, ur_longitude, ur_lat, ur_lon;
324
325
	GRIB_CHECK(grib_get_double(h,"latitudeOfLastGridPointInDegrees",&ur_latitude),0);
	GRIB_CHECK(grib_get_double(h,"longitudeOfLastGridPointInDegrees",&ur_longitude),0);
326
	Coords::rotatedToTrueLatLon(latitudeOfNorthernPole, longitudeOfNorthernPole, ur_latitude, ur_longitude, ur_lat, ur_lon);
327
	double cntr_lat, cntr_lon; //geographic coordinates
328
	Coords::rotatedToTrueLatLon(latitudeOfNorthernPole, longitudeOfNorthernPole, .5*(ll_latitude+ur_latitude), .5*(ll_longitude+ur_longitude), cntr_lat, cntr_lon);
329
330

	//determining cell size
331
332
333
	long Ni, Nj;
	GRIB_CHECK(grib_get_long(h,"Ni",&Ni),0);
	GRIB_CHECK(grib_get_long(h,"Nj",&Nj),0);
334
	double bearing;
335
336
337
338
339
340
341
342
	cell_x = Coords::VincentyDistance(cntr_lat, ll_lon, cntr_lat, ur_lon, bearing) / (double)Ni;
	cell_y = Coords::VincentyDistance(cntr_lon, ll_lat, cntr_lon, ur_lat, bearing) / (double)Nj;

	//determining bearing offset
	double delta_lat, delta_lon; //geographic coordinates
	Coords::rotatedToTrueLatLon(latitudeOfNorthernPole, longitudeOfNorthernPole, .5*(ll_latitude+ur_latitude)+1., .5*(ll_longitude+ur_longitude), delta_lat, delta_lon);
	Coords::VincentyDistance(cntr_lat, cntr_lon, delta_lat, delta_lon, bearing_offset);
	bearing_offset = fmod( bearing_offset + 180., 360.) - 180.; // turn into [-180;180)
343

344
345
346
	//computing lower left corner by using the center point as reference
	Coords cntr(coordin, coordinparam);
	cntr.setLatLon(cntr_lat, cntr_lon, IOUtils::nodata);
347
348
	const double cellsize=.5*(cell_x+cell_y);
	cntr.moveByXY(-.5*(double)Ni*cellsize, -.5*(double)Nj*cellsize);
349

350
	//checking that cellsize does not vary too much across the grid
351
	/*const double cellsize_x_ll = Coords::lon_degree_lenght(ll_latitude)*d_j;
352
353
354
355
356
357
	const double cellsize_x_ur = Coords::lon_degree_lenght(ur_latitude)*d_j;
	if( fabs(cellsize_x_ll-cellsize_x_ur)/cellsize_x > 1./100.) {
		stringstream ss;
		ss << "Cell size varying too much in the x direction between lower left and upper right corner: ";
		ss << cellsize_x_ll << "m to " << cellsize_x_ur << "m";
		throw IOException(ss.str(), AT);
358
	}*/
359

360
	return cntr; //this is now the ll corner
361
362
}

363
void GRIBIO::read2Dlevel(grib_handle* h, Grid2DObject& grid_out, const bool& read_geolocalization)
364
{
365
366
367
368
369
370
371
372
	long Ni, Nj;
	GRIB_CHECK(grib_get_long(h,"Ni",&Ni),0);
	GRIB_CHECK(grib_get_long(h,"Nj",&Nj),0);

	double *values;
	size_t values_len= 0;
	GRIB_CHECK(grib_get_size(h,"values",&values_len),0);
	if(values_len!=(unsigned)(Ni*Nj)) {
373
		stringstream ss;
374
375
376
377
378
379
380
		ss << "Declaring grid of size " << Ni << "x" << Nj << "=" << Ni*Nj << " ";
		ss << "but containing " << values_len << " values. This is inconsistent!";
		throw InvalidArgumentException(ss.str(), AT);
	}
	values = (double*)malloc(values_len*sizeof(double));

	GRIB_CHECK(grib_get_double_array(h,"values",values,&values_len),0);
381
382
383
384
385
386

	if(read_geolocalization) {
		llcorner = getGeolocalization(h, cellsize_x, cellsize_y);
		if( fabs(cellsize_x-cellsize_y)/cellsize_x > 1./100.) {
			throw InvalidArgumentException("Cells can not be represented by square cells. This is not supported!", AT);
		}
387
	}
388
	grid_out.set(static_cast<unsigned int>(Ni), static_cast<unsigned int>(Nj), .5*(cellsize_x+cellsize_y), llcorner);
389
390
391
392
	int i=0;
	for(unsigned int jj=0; jj<(unsigned)Nj; jj++) {
		for(unsigned int ii=0; ii<(unsigned)Ni; ii++)
			grid_out(ii,jj) = values[i++];
393
	}
394
395
	free(values);
}
396

397
bool GRIBIO::read2DGrid_indexed(const double& in_marsParam, const long& i_levelType, const long& i_level, const Date i_date, Grid2DObject& grid_out)
398
{
399
400
401
	GRIB_CHECK(grib_index_select_double(idx,"marsParam",in_marsParam),0);
	GRIB_CHECK(grib_index_select_long(idx,"indicatorOfTypeOfLevel", i_levelType),0);

402
403
	grib_handle* h=NULL;
	int err=0;
404
	while((h = grib_handle_new_from_index(idx,&err)) != NULL) {
405
406
		if(!h) {
			cleanup();
407
			throw IOException("Unable to create grib handle from index", AT);
408
409
		}

410
411
412
413
414
415
416
417
		Date base_date;
		double P1, P2;
		getDate(h, base_date, P1, P2);

		//see WMO code table5 for definitions of timeRangeIndicator. http://dss.ucar.edu/docs/formats/grib/gribdoc/timer.html
		long timeRange;
		GRIB_CHECK(grib_get_long(h,"timeRangeIndicator", &timeRange),0);

418
		long level=0;
419
		if(i_level!=0) GRIB_CHECK(grib_get_long(h,"level", &level),0);
420
421
422
423
424
425
426
427
428
429
		if(level==i_level) {
			//geolocalization has been initialized when indexing the file, so we don't need to redo it
			if( (i_date.isUndef()) ||
			    (timeRange==0 && i_date==base_date+P1) ||
			    (timeRange==1 && i_date==base_date) ||
			    ((timeRange==2 || timeRange==3) && i_date>=base_date+P1 && i_date<=base_date+P2) ||
			    ((timeRange==4 || timeRange==5) && i_date==base_date+P2) ) {
				read2Dlevel(h, grid_out, false);
				return true;
			}
430
		}
431
		grib_handle_delete(h);
432
	}
433
	return false;
434
435
}

436
void GRIBIO::read2DGrid(Grid2DObject& grid_out, const std::string& i_name)
437
{
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
	const std::string filename = grid2dpath_in+"/"+i_name;
	fp = fopen(filename.c_str(),"r");
	if(fp==NULL) {
		stringstream ss;
		ss << "Error openning file \"" << filename << "\", possible reason: " << strerror(errno);
		throw FileAccessException(ss.str(), AT);
	}

	grib_handle* h=NULL;
	int err=0;
	if((h = grib_handle_new_from_file(0,fp,&err)) != NULL) {
		if(!h) {
			cleanup();
			throw IOException("Unable to create grib handle for \""+filename+"\"", AT);
		}

454
		read2Dlevel(h, grid_out, true);
455
456
457
458
459
460
461
		grib_handle_delete(h);
	} else {
		cleanup();
		throw IOException("No grid found in file \""+filename+"\"", AT);
	}

	cleanup();
462
463
}

464
void GRIBIO::indexFile(const std::string& filename)
465
{
466
467
468
469
470
471
472
	fp = fopen(filename.c_str(),"r");
	if(fp==NULL) {
		stringstream ss;
		ss << "Error openning file \"" << filename << "\", possible reason: " << strerror(errno);
		throw FileAccessException(ss.str(), AT);
	}

473
	int err=0;
474
	std::string keys("marsParam:d,indicatorOfTypeOfLevel:l"); //indexing keys
475
476
477
478
	char *c_filename = (char *)filename.c_str();
	idx = grib_index_new_from_file(0, c_filename, keys.c_str(), &err);
	if(err!=0) {
		cleanup();
479
		throw IOException("Failed to index GRIB file \""+filename+"\". Is it a valid GRIB file?", AT);
480
481
482
	}
	indexed=true;
	idx_filename=filename;
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502

	//read geolocalization of the first grid we find
	grib_handle* h=NULL;
	err=0;
	if((h = grib_handle_new_from_file(0,fp,&err)) != NULL) {
		if(!h) {
			cleanup();
			throw IOException("Unable to create grib handle for \""+filename+"\"", AT);
		}

		llcorner = getGeolocalization(h, cellsize_x, cellsize_y); //this sets llcorner, cellsize and bearing_offset
		if( fabs(cellsize_x-cellsize_y)/cellsize_x > 1./100.) {
			throw InvalidArgumentException("Cells can not be represented by square cells. This is not supported!", AT);
		}
		grib_handle_delete(h);
	} else {
		cleanup();
		throw IOException("No grid found in file \""+filename+"\"", AT);
	}

503
}
504

505
void GRIBIO::read2DGrid(Grid2DObject& grid_out, const MeteoGrids::Parameters& parameter, const Date& date)
506
{
507
508
	Date UTC_date = date;
	UTC_date.setTimeZone(tz_in);
509

510
	const std::string filename = grid2dpath_in+"/"+prefix+UTC_date.toString(Date::NUM).substr(0,10)+ext;
511
512

	read2DGrid(filename, grid_out, parameter, UTC_date);
513
514
}

515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
void GRIBIO::readWind(const std::string& filename, const Date& date)
{
	if(wind_date==date) return; //wind fields are already up to date

	if(read2DGrid_indexed(32.2, 105, 10, date, VW)) { //FF_10M
		if(!read2DGrid_indexed(31.2, 105, 10, date, DW)) //DD_10M
			throw NoAvailableDataException("Can not read wind direction in file \""+filename+"\"", AT);
	} else {
		Grid2DObject U,V;
		read2DGrid_indexed(33.2, 105, 10, date, U); //U_10M, also in 110, 10 as U
		read2DGrid_indexed(34.2, 105, 10, date, V); //V_10M, also in 110, 10 as V

		VW.set(U.ncols, U.nrows, U.cellsize, U.llcorner);
		DW.set(U.ncols, U.nrows, U.cellsize, U.llcorner);
		for(unsigned int jj=0; jj<VW.nrows; jj++) {
			for(unsigned int ii=0; ii<VW.ncols; ii++) {
				VW(ii,jj) = sqrt( IOUtils::pow2(U(ii,jj)) + IOUtils::pow2(V(ii,jj)) );
				DW(ii,jj) = fmod( atan2( U(ii,jj), V(ii,jj) ) * to_deg + 360. + bearing_offset, 360.); // turn into degrees [0;360)
			}
		}
	}

	wind_date = date;
}

540
void GRIBIO::read2DGrid(const std::string& filename, Grid2DObject& grid_out, const MeteoGrids::Parameters& parameter, const Date& date)
541
{ //Parameters should be read in table 2 if available since this table is the one standardized by WMO
542
	if(!indexed || idx_filename!=filename) {
543
544
545
546
		cleanup();
		indexFile(filename);
	}

547
	//Basic meteo parameters
548
	if(parameter==MeteoGrids::P) read2DGrid_indexed(1.2, 1, 0, date, grid_out); //PS
549
	if(parameter==MeteoGrids::TA) read2DGrid_indexed(11.2, 105, 2, date, grid_out); //T_2M
550
	if(parameter==MeteoGrids::RH) {
551
		if(!read2DGrid_indexed(52.2, 105, 2, date, grid_out)) { //RELHUM_2M
552
			Grid2DObject ta;
553
554
			read2DGrid_indexed(11.2, 105, 2, date, ta); //T_2M
			read2DGrid_indexed(17.2, 105, 2, date, grid_out); //TD_2M
555
556
557
558
559
560
561
			for(unsigned int jj=0; jj<grid_out.nrows; jj++) {
				for(unsigned int ii=0; ii<grid_out.ncols; ii++) {
					grid_out(ii,jj) = Atmosphere::DewPointtoRh(grid_out(ii,jj), ta(ii,jj), true);
				}
			}
		}
	}
562
563
	if(parameter==MeteoGrids::TSS) read2DGrid_indexed(197.201, 111, 0, date, grid_out); //T_SO
	if(parameter==MeteoGrids::TSG) read2DGrid_indexed(11.2, 1, 0, date, grid_out); //T_G
564
565

	//hydrological parameters
566
567
568
569
	if(parameter==MeteoGrids::HNW) read2DGrid_indexed(61.2, 1, 0, date, grid_out); //tp
	if(parameter==MeteoGrids::ROT) read2DGrid_indexed(90.2, 112, 0, date, grid_out); //RUNOFF
	if(parameter==MeteoGrids::SWE) read2DGrid_indexed(65.2, 1, 0, date, grid_out); //W_SNOW
	if(parameter==MeteoGrids::HS) {
570
571
572
573
574
575
		if(!read2DGrid_indexed(66.2, 1, 0, date, grid_out)) {
			Grid2DObject snow_density;
			read2DGrid_indexed(133.201, 1, 0, date, snow_density); //RHO_SNOW
			read2DGrid_indexed(65.2, 1, 0, date, grid_out); //W_SNOW
			grid_out.grid2D /= snow_density.grid2D;
		}
576
	}
577
578

	//radiation parameters
579
580
581
582
	if(parameter==MeteoGrids::ALB) {
		read2DGrid_indexed(84.2, 1, 0, date, grid_out); //ALB_RAD
		grid_out.grid2D /= 100.;
	}
583
584
585
586
587
	if(parameter==MeteoGrids::ILWR) {
		if(read2DGrid_indexed(115.2, 1, 0, date, grid_out)) { //long wave
			grid_out.grid2D *= -1.;
		} else read2DGrid_indexed(25.201, 1, 0, date, grid_out); //ALWD_S
	}
588
	/*if(parameter==MeteoGrids::CLD) { //cloudiness
589
590
591
592
		if(read2DGrid_indexed(74.2, 1, 0, date, grid_out)) //CLCM
		grid_out.grid2D /= 100.;
	}*/

593
	if(parameter==MeteoGrids::ISWR) {
594
595
596
597
		if(read2DGrid_indexed(116.2, 1, 0, date, grid_out)) { //short wave
			grid_out.grid2D *= -1.;
		} else {
		//if(!read2DGrid_indexed(111.250, 1, 0, date, grid_out)) { //GLOB
598
			Grid2DObject diff;
599
600
			read2DGrid_indexed(23.201, 1, 0, date, diff); //diffuse rad, ASWDIFD_S
			read2DGrid_indexed(22.201, 1, 0, date, grid_out); //direct rad, ASWDIR_S
601
602
603
604
605
			grid_out.grid2D += diff.grid2D;
		}
	}

	//DEM parameters
606
	if(parameter==MeteoGrids::DEM) read2DGrid_indexed(8.2, 1, 0, date, grid_out); //HSURF
607
608
609
610
611
612
	if(parameter==MeteoGrids::SLOPE) {
		read2DGrid_indexed(98.202, 1, 0, date, grid_out); //SLO_ANG
		grid_out.grid2D *= to_deg;
	}
	if(parameter==MeteoGrids::AZI) {
		read2DGrid_indexed(99.202, 1, 0, date, grid_out); //SLO_ASP
613
614
615
616
617
		for(unsigned int jj=0; jj<grid_out.nrows; jj++) {
			for(unsigned int ii=0; ii<grid_out.ncols; ii++) {
				grid_out(ii,jj) = fmod( grid_out(ii,jj)*to_deg + 360. + bearing_offset, 360.); // turn into degrees [0;360)
			}
		}
618
	}
619

620
	//Wind parameters
621
	if(parameter==MeteoGrids::VW_MAX) read2DGrid_indexed(187.201, 105, 10, date, grid_out); //VMAX_10M 10m
622
623
624
625
626
627
628
	if(parameter==MeteoGrids::W) read2DGrid_indexed(40.2, 109, 10, date, grid_out); //W, 10m
	 //we need to use VW, DW, correct for re-projection and recompute U,V
	if(parameter==MeteoGrids::U) {
		readWind(filename, date);
		for(unsigned int jj=0; jj<grid_out.nrows; jj++) {
			for(unsigned int ii=0; ii<grid_out.ncols; ii++) {
				grid_out(ii,jj) = VW(ii,jj)*sin(DW(ii,jj)*to_rad);
629
630
631
			}
		}
	}
632
633
634
635
636
	if(parameter==MeteoGrids::V) {
		readWind(filename, date);
		for(unsigned int jj=0; jj<grid_out.nrows; jj++) {
			for(unsigned int ii=0; ii<grid_out.ncols; ii++) {
				grid_out(ii,jj) = VW(ii,jj)*cos(DW(ii,jj)*to_rad);
637
638
639
			}
		}
	}
640
641
642
643
644
645
646
647
	if(parameter==MeteoGrids::DW) {
		readWind(filename, date);
		grid_out = DW;
	}
	if(parameter==MeteoGrids::VW) {
		readWind(filename, date);
		grid_out = VW;
	}
648
649
650
651
652
653
654
655
656

	if(grid_out.isEmpty()) {
		stringstream ss;
		ss << "No suitable data found for parameter " << MeteoGrids::getParameterName(parameter) << " ";
		ss << "at time step " << date.toString(Date::ISO) << " in file \"" << filename << "\"";
		throw NoAvailableDataException(ss.str(), AT);
	}

	//correcting wind speeds
657
658
	/*if(parameter==MeteoGrids::U || parameter==MeteoGrids::V || parameter==MeteoGrids::W || parameter==MeteoGrids::VW || parameter==MeteoGrids::VW_MAX) {
		//we need to compute the wind at 7.5m
659
		Grid2DObject Z0;
660
		if(read2DGrid_indexed(83.2, 1, 0, date, Z0)) { //Z0
661
662
			for(unsigned int jj=0; jj<grid_out.nrows; jj++) {
				for(unsigned int ii=0; ii<grid_out.ncols; ii++) {
663
					grid_out(ii,jj) = Atmosphere::windLogProfile(grid_out(ii,jj), 10., 7.5, Z0(ii,jj));
664
665
				}
			}
666
		} else {
667
			const double wind_factor = Atmosphere::windLogProfile(1., 10., 7.5, 0.03);
668
			grid_out.grid2D *= wind_factor;
669
		}
670
	}*/
671
672
}

673
void GRIBIO::readDEM(DEMObject& dem_out)
674
{
675
	const Date d; //ie: undef. This will be caught when reading the GRIB file
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
	std::string filename;

	cfg.getValue("DEMFILE", "Input", filename);
	read2DGrid(filename, dem_out, MeteoGrids::DEM, d);
	if(update_dem) {
		dem_out.update();
	} else {
		const int dem_ppt=dem_out.getUpdatePpt();
		if(dem_ppt&DEMObject::SLOPE) {
			Grid2DObject slope;
			read2DGrid(filename, slope, MeteoGrids::SLOPE, d);
			dem_out.slope=slope.grid2D;
			Grid2DObject azi;
			read2DGrid(filename, azi, MeteoGrids::AZI, d);
			dem_out.azi=azi.grid2D;
		}
		if(dem_ppt&DEMObject::NORMAL || dem_ppt&DEMObject::CURVATURE) {
			//we will only update the normals and/or curvatures, then revert update properties
			if(dem_ppt&DEMObject::NORMAL && dem_ppt&DEMObject::CURVATURE) dem_out.setUpdatePpt((DEMObject::update_type)(DEMObject::NORMAL|DEMObject::CURVATURE));
			else if(dem_ppt&DEMObject::NORMAL) dem_out.setUpdatePpt(DEMObject::NORMAL);
			else if(dem_ppt&DEMObject::CURVATURE) dem_out.setUpdatePpt(DEMObject::CURVATURE);

			dem_out.update();
			dem_out.setUpdatePpt((DEMObject::update_type)dem_ppt);
		}

		dem_out.updateAllMinMax();
	}
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
}

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

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

void GRIBIO::readStationData(const Date&, std::vector<StationData>& /*vecStation*/)
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

724
725
726
727
728
729
730
731
732
733
734
735
void GRIBIO::scanMeteoPath()
{
	std::list<std::string> dirlist;
	IOUtils::readDirectory(meteopath_in, dirlist, ext);
	dirlist.sort();

	//Check date in every filename and cache it
	std::list<std::string>::iterator it = dirlist.begin();
	while ((it != dirlist.end())) {
		const std::string& filename = *it;
		std::string::size_type spos = filename.find_first_of("0123456789");
		Date date;
736
		IOUtils::convertString(date, filename.substr(spos,10), tz_in);
737
738
739
740
741
742
743
		std::pair<Date,std::string> tmp(date, filename);

		cache_meteo_files.push_back(tmp);
		it++;
	}
}

744
745
void GRIBIO::readMeteoData(const Date& dateStart, const Date& dateEnd,
                             std::vector< std::vector<MeteoData> >& vecMeteo,
746
747
                             const size_t&)
{
748
749
750
751
752
	if(!meteo_initialized) {
		readStations();
		scanMeteoPath();
		meteo_initialized=true;
	}
753

754
	vecMeteo.clear();
755

756
757
	double *lats = (double*)malloc(vecPts.size()*sizeof(double));
	double *lons = (double*)malloc(vecPts.size()*sizeof(double));
758
759
760
761
762
763
764
765
766
767
768
769
770
771
	std::vector<StationData> meta; //metadata for meteo time series
	bool meta_ok=false; //set to true once the metadata have been read

	//find index of first time step
	unsigned int idx_start;
	bool start_found=false;
	for(idx_start=0; idx_start<cache_meteo_files.size(); idx_start++) {
		if(dateStart<cache_meteo_files[idx_start].first) {
			start_found=true;
			break;
		}
	}
	if(start_found==false) return;
	if(idx_start>0) idx_start--; //start with first element before dateStart (useful for resampling)
772
773

	try {
774
775
776
777
		for(unsigned int ii=idx_start; ii<cache_meteo_files.size(); ii++) {
			const Date& date = cache_meteo_files[ii].first;
			if(date>dateEnd) break;
			const std::string filename = meteopath_in+"/"+cache_meteo_files[ii].second;
778
779
780

			if(!indexed || idx_filename!=filename) {
				cleanup();
781
				indexFile(filename); //this will also read geolocalization
782
783
			}
			if(!meta_ok) {
784
				if(readMeteoMeta(vecPts, meta, lats, lons)==false) {
785
786
787
788
					//some points have been removed vecPts has been changed -> re-reading
					free(lats); free(lons);
					lats = (double*)malloc(vecPts.size()*sizeof(double));
					lons = (double*)malloc(vecPts.size()*sizeof(double));
789
					readMeteoMeta(vecPts, meta, lats, lons);
790
				}
791
				vecMeteo.insert(vecMeteo.begin(), vecPts.size(), std::vector<MeteoData>()); //allocation for the vectors now that we know how many true stations we have
792
				meta_ok=true;
793
794
795
796
797
798
799
800
801
			}

			std::vector<MeteoData> Meteo;
			readMeteoStep(meta, lats, lons, date, Meteo);
			for(unsigned int jj=0; jj<vecPts.size(); jj++)
				vecMeteo[jj].push_back(Meteo[jj]);
		}
	} catch(...) {
		free(lats); free(lons);
802
		cleanup();
803
		throw;
804
805
	}

806
	free(lats); free(lons);
807
808
}

809
810
811
812
813
814
815
816
817
bool GRIBIO::removeDuplicatePoints(std::vector<Coords>& vecPts, double *lats, double *lons)
{ //remove potential duplicates. Returns true if some have been removed
	const unsigned int npoints = vecPts.size();
	std::vector<size_t> deletions;
	deletions.reserve(npoints);
	for(unsigned int ii=0; ii<npoints; ii++) {
		const double lat = lats[ii];
		const double lon = lons[ii];
		for(unsigned int jj=ii+1; jj<npoints; jj++) {
818
819
820
			if(lat==lats[jj] && lon==lons[jj]) {
				deletions.push_back(jj);
			}
821
822
823
824
825
826
		}
	}

	//we need to erase from the end in order to keep the index unchanged...
	for(unsigned int ii=deletions.size(); ii>0; ii--) {
		const unsigned int index=deletions[ii-1];
827
		vecPts.erase(vecPts.begin()+index);
828
829
830
831
832
833
	}

	if(deletions.size()>0) return true;
	return false;
}

834
bool GRIBIO::readMeteoMeta(std::vector<Coords>& vecPts, std::vector<StationData> &stations, double *lats, double *lons)
835
836
837
{//return true if the metadata have been read, false if it needs to be re-read (ie: some points were leading to duplicates -> vecPts has been changed)
	stations.clear();

838
839
840
841
	GRIB_CHECK(grib_index_select_double(idx,"marsParam",8.2),0); //This is the DEM
	GRIB_CHECK(grib_index_select_long(idx,"indicatorOfTypeOfLevel", 1),0);

	int err=0;
842
843
844
845
846
	grib_handle* h = grib_handle_new_from_index(idx,&err);
	if(h==NULL) {
		cleanup();
		throw IOException("Can not find DEM grid in GRIB file!", AT);
	}
847

848
	const long npoints = vecPts.size();
849
	double latitudeOfSouthernPole, longitudeOfSouthernPole;
850
851
	GRIB_CHECK(grib_get_double(h,"latitudeOfSouthernPoleInDegrees",&latitudeOfSouthernPole),0);
	GRIB_CHECK(grib_get_double(h,"longitudeOfSouthernPoleInDegrees",&longitudeOfSouthernPole),0);
852
853
	latitudeOfNorthernPole = -latitudeOfSouthernPole;
	longitudeOfNorthernPole = longitudeOfSouthernPole+180.;
854

855
856
	long Ni;
	GRIB_CHECK(grib_get_long(h,"Ni",&Ni),0);
857

858
859
	//build GRIB local coordinates for the points
	for(unsigned int ii=0; ii<(unsigned)npoints; ii++) {
860
		Coords::trueLatLonToRotated(latitudeOfNorthernPole, longitudeOfNorthernPole, vecPts[ii].getLat(), vecPts[ii].getLon(), lats[ii], lons[ii]);
861
	}
862

863
864
865
866
867
868
869
870
871
872
873
	//retrieve nearest points
	double *outlats = (double*)malloc(npoints*sizeof(double));
	double *outlons = (double*)malloc(npoints*sizeof(double));
	double *values = (double*)malloc(npoints*sizeof(double));
	double *distances = (double*)malloc(npoints*sizeof(double));
	int *indexes = (int *)malloc(npoints*sizeof(int));
	if(grib_nearest_find_multiple(h, 0, lats, lons, npoints, outlats, outlons, values, distances, indexes)!=0) {
		grib_handle_delete(h);
		cleanup();
		throw IOException("Errro when searching for nearest points in DEM", AT);
	}
874

875
876
	//remove potential duplicates
	if(removeDuplicatePoints(vecPts, outlats, outlons)==true) {
877
878
		free(outlats); free(outlons); free(values); free(distances); free(indexes);
		grib_handle_delete(h);
879
		return false;
880
	}
881
882
883
884
885

	//fill metadata
	for(unsigned int ii=0; ii<(unsigned)npoints; ii++) {
		StationData sd;
		sd.position.setProj(coordin, coordinparam);
886
		double true_lat, true_lon;
887
		Coords::rotatedToTrueLatLon(latitudeOfNorthernPole, longitudeOfNorthernPole, outlats[ii], outlons[ii], true_lat, true_lon);
888
		sd.position.setLatLon(true_lat, true_lon, values[ii]);
889
890
891
892
893
894
895
896
897
898
899
900
		stringstream ss;
		ss << "Point_" << indexes[ii];
		sd.stationID=ss.str();
		stringstream ss2;
		ss2 << "GRIB point (" << indexes[ii] % Ni << "," << indexes[ii] / Ni << ")";
		sd.stationName=ss2.str();
		stations.push_back(sd);
	}

	free(outlats); free(outlons); free(values); free(distances); free(indexes);
	grib_handle_delete(h);
	return true;
901
902
903
904
905
906
}

bool GRIBIO::readMeteoValues(const double& marsParam, const long& levelType, const long& i_level, const Date& i_date, const long& npoints, double *lats, double *lons, double *values)
{
	GRIB_CHECK(grib_index_select_double(idx,"marsParam",marsParam),0);
	GRIB_CHECK(grib_index_select_long(idx,"indicatorOfTypeOfLevel", levelType),0);
907
908
909
910
911
912
913
914
915

	grib_handle* h=NULL;
	int err=0;
	while((h = grib_handle_new_from_index(idx,&err)) != NULL) {
		if(!h) {
			cleanup();
			throw IOException("Unable to create grib handle from index", AT);
		}

916
917
918
919
920
921
922
923
		Date base_date;
		double P1, P2;
		getDate(h, base_date, P1, P2);

		//see WMO code table5 for definitions of timeRangeIndicator. http://dss.ucar.edu/docs/formats/grib/gribdoc/timer.html
		long timeRange;
		GRIB_CHECK(grib_get_long(h,"timeRangeIndicator", &timeRange),0);

924
925
		long level=0;
		if(i_level!=0) GRIB_CHECK(grib_get_long(h,"level", &level),0);
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
		if(level==i_level) {
			if( (i_date.isUndef()) ||
			    (timeRange==0 && i_date==base_date+P1) ||
			    (timeRange==1 && i_date==base_date) ||
			    ((timeRange==2 || timeRange==3) && i_date>=base_date+P1 && i_date<=base_date+P2) ||
			    ((timeRange==4 || timeRange==5) && i_date==base_date+P2) ) {
				double *outlats = (double*)malloc(npoints*sizeof(double));
				double *outlons = (double*)malloc(npoints*sizeof(double));
				double *distances = (double*)malloc(npoints*sizeof(double));
				int *indexes = (int *)malloc(npoints*sizeof(int));
				if(grib_nearest_find_multiple(h, 0, lats, lons, npoints, outlats, outlons, values, distances, indexes)!=0) {
					grib_handle_delete(h);
					cleanup();
					throw IOException("Errro when searching for nearest points in DEM", AT);
				}

				free(outlats); free(outlons); free(distances); free(indexes);
943
				grib_handle_delete(h);
944
				return true;
945
946
947
948
949
			}
		}
		grib_handle_delete(h);
	}
	return false;
950
951
}

952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
void GRIBIO::fillMeteo(double *values, const MeteoData::Parameters& param, const long& npoints, std::vector<MeteoData> &Meteo) {
	for(unsigned int ii=0; ii<(unsigned)npoints; ii++) {
		Meteo[ii](param) = values[ii];
	}
}

void GRIBIO::readMeteoStep(std::vector<StationData> &stations, double *lats, double *lons, const Date i_date, std::vector<MeteoData> &Meteo)
{
	const long npoints = stations.size();

	for(unsigned int ii=0; ii<(unsigned)npoints; ii++) {
		MeteoData md;
		md.meta = stations[ii];
		md.date = i_date;
		Meteo.push_back(md);
	}

	double *values = (double*)malloc(npoints*sizeof(double));
	double *values2 = (double*)malloc(npoints*sizeof(double)); //for extra parameters

	//basic meteorological parameters
	if(readMeteoValues(1.2, 1, 0, i_date, npoints, lats, lons, values)) fillMeteo(values, MeteoData::P, npoints, Meteo); //PS
	if(readMeteoValues(11.2, 105, 2, i_date, npoints, lats, lons, values)) fillMeteo(values, MeteoData::TA, npoints, Meteo); //T_2M
	if(readMeteoValues(197.201, 111, 0, i_date, npoints, lats, lons, values)) fillMeteo(values, MeteoData::TSS, npoints, Meteo); //T_SO
	if(readMeteoValues(11.2, 1, 0, i_date, npoints, lats, lons, values)) fillMeteo(values, MeteoData::TSG, npoints, Meteo); //T_G
	if(readMeteoValues(52.2, 105, 2, i_date, npoints, lats, lons, values)) fillMeteo(values, MeteoData::RH, npoints, Meteo); //RELHUM_2M
	else if(readMeteoValues(17.2, 105, 2, i_date, npoints, lats, lons, values)) { //TD_2M
		for(unsigned int ii=0; ii<(unsigned)npoints; ii++) {
			if(Meteo[ii](MeteoData::TA)!=IOUtils::nodata)
				Meteo[ii](MeteoData::RH) = Atmosphere::DewPointtoRh(values[ii], Meteo[ii](MeteoData::TA), true);
		}
	}

	//hydrological parameters
	if(readMeteoValues(61.2, 1, 0, i_date, npoints, lats, lons, values)) fillMeteo(values, MeteoData::HNW, npoints, Meteo); //tp
987
988
	if(readMeteoValues(66.2, 1, 0, i_date, npoints, lats, lons, values)) fillMeteo(values, MeteoData::HS, npoints, Meteo);
	else if(readMeteoValues(133.201, 1, 0, i_date, npoints, lats, lons, values)  //RHO_SNOW
989
990
991
992
993
994
995
	   && readMeteoValues(65.2, 1, 0, i_date, npoints, lats, lons, values2)) { //W_SNOW
		for(unsigned int ii=0; ii<(unsigned)npoints; ii++) {
			Meteo[ii](MeteoData::HS) = values2[ii] / values[ii];
		}
	}

	//radiation parameters
996
997
998
999
1000
	if(readMeteoValues(115.2, 1, 0, i_date, npoints, lats, lons, values)) { //long wave
		for(unsigned int ii=0; ii<(unsigned)npoints; ii++) {
			Meteo[ii](MeteoData::ISWR) = -values[ii];
		}
	} else if(readMeteoValues(25.201, 1, 0, i_date, npoints, lats, lons, values)) fillMeteo(values, MeteoData::ILWR, npoints, Meteo); //ALWD_S