WSL/SLF GitLab Repository

GRIBIO.cc 41.5 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
#include <meteoio/MathOptim.h>
24

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

using namespace std;

namespace mio {
/**
33
 * @page gribio GRIBIO
34
 * @section gribio_format Format and limitations
35
 * This plugin reads GRIB (https://en.wikipedia.org/wiki/GRIB) files as produced by meteorological models.
36
 * 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)).
37
 * 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
38
 * (levels description is available at http://www.nco.ncep.noaa.gov/pmb/docs/on388/). Standard COSMO grids are listed at http://zephyr.ucd.ie/mediawiki/index.php/COSMO_GRIB .
39
 *
40
41
42
43
44
45
46
47
48
 * 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.
49
 *
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
 * @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
 *
65
66
 * @section gribio_units Units
 * As specified by WMO.
67
 *
68
 * @section gribio_keywords Keywords
69
 * This plugin uses the following keywords:
70
71
72
 * - COORDSYS: coordinate system (see Coords)
 * - COORDPARAM: extra coordinates parameters (see Coords)
 * - GRID2DPATH: path where to find the grids
73
74
 * - GRID2DPREFIX: prefix to append when generating a file name for reading (ie: something like "laf" for Cosmo-Analysis-full domain), optional
 * - GRID2DEXT: grib file extension, or <i>none</i> for no file extension (default: .grb)
75
76
 * - 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)
77
78
 * - METEOPATH: path where to find the grids for extracting time series at special points
 * - METEOEXT: file extension, or <i>none</i> for no file extension (default: .grb)
79
 * - STATION#: coordinates for virtual stations (if using GRIB as METEO plugin). Each station is given by its coordinates and the closest
80
81
 * 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.
82
 *
83
84
85
 */

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

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

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

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

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

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

	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);
	}
142
	cfg.getValue("GRID2DPREFIX", "Input", grid2d_prefix, Config::nothrow);
143

144
145
146
147
148
149
150
	meteo_ext = default_ext;
	cfg.getValue("METEOEXT", "Input", meteo_ext, Config::nothrow);
	if(meteo_ext=="none") meteo_ext="";

	grid2d_ext = default_ext;
	cfg.getValue("GRID2DEXT", "Input", grid2d_ext, Config::nothrow);
	if(grid2d_ext=="none") grid2d_ext="";
151
}
152

153
154
155
156
157
158
159
160
161
162
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);
163
		IOUtils::stripComments(current_station);
164
165
166

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

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);
198
199
}

200
void GRIBIO::listKeys(grib_handle** h, const std::string& filename)
201
{
202
	const unsigned int MAX_VAL_LEN=1024; //max value string length in GRIB
203
	unsigned long key_iterator_filter_flags=GRIB_KEYS_ITERATOR_ALL_KEYS;
204
205
206
207
208
209
210
	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);
	}
211

212
213
214
215
216
217
218
219
	//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);
220
		std::cerr << name << " = " << value << "\n";
221
222
223
224
225
	}

	grib_keys_iterator_delete(kiter);
}

226
void  GRIBIO::listFields(const std::string& filename)
227
{
228
	grib_handle* h=NULL;
229
	int err=0;
230

231
	//loop over the messages
232
233
234
235
236
237
	while((h = grib_handle_new_from_file(0,fp,&err)) != NULL) {
		if(!h) {
			cleanup();
			throw IOException("Unable to create grib handle for \""+filename+"\"", AT);
		}

238
239
240
241
242
243
244
245
		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);
246
		long levelType;
247
248
249
		GRIB_CHECK(grib_get_long(h,"indicatorOfTypeOfLevel", &levelType),0);
		long level=0;
		if(levelType!=1) GRIB_CHECK(grib_get_long(h,"level", &level),0);
250
		std::cerr << marsParam << " " << shortname << " " << name << " type " << levelType << " level " << level << "\n";
251
		grib_handle_delete(h);
252
253
	}
}
254

255
void GRIBIO::getDate(grib_handle* h, Date &base, double &d1, double &d2) {
256
257
258
259
260
261
	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;
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
302
303
304
305
306
	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;
307
308
}

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

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

	//determining cell size
336
337
338
	long Ni, Nj;
	GRIB_CHECK(grib_get_long(h,"Ni",&Ni),0);
	GRIB_CHECK(grib_get_long(h,"Nj",&Nj),0);
339
	double bearing;
340
341
342
343
344
345
346
347
	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)
348

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

355
	//checking that cellsize does not vary too much across the grid
356
	/*const double cellsize_x_ll = Coords::lon_degree_lenght(ll_latitude)*d_j;
357
358
359
360
361
362
	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);
363
	}*/
364

365
	return cntr; //this is now the ll corner
366
367
}

368
void GRIBIO::read2Dlevel(grib_handle* h, Grid2DObject& grid_out, const bool& read_geolocalization)
369
{
370
371
372
373
374
375
376
377
	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)) {
378
		stringstream ss;
379
380
381
382
383
384
385
		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);
386
387
388
389
390
391

	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);
		}
392
	}
393
	grid_out.set(static_cast<unsigned int>(Ni), static_cast<unsigned int>(Nj), .5*(cellsize_x+cellsize_y), llcorner);
394
395
396
397
	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++];
398
	}
399
400
	free(values);
}
401

402
bool GRIBIO::read2DGrid_indexed(const double& in_marsParam, const long& i_levelType, const long& i_level, const Date i_date, Grid2DObject& grid_out)
403
{
404
405
406
	GRIB_CHECK(grib_index_select_double(idx,"marsParam",in_marsParam),0);
	GRIB_CHECK(grib_index_select_long(idx,"indicatorOfTypeOfLevel", i_levelType),0);

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

415
416
417
418
419
420
421
422
		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);

423
		long level=0;
424
		if(i_level!=0) GRIB_CHECK(grib_get_long(h,"level", &level),0);
425
426
427
428
429
430
431
432
433
434
		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;
			}
435
		}
436
		grib_handle_delete(h);
437
	}
438
	return false;
439
440
}

441
void GRIBIO::read2DGrid(Grid2DObject& grid_out, const std::string& i_name)
442
{
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
	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);
		}

459
		read2Dlevel(h, grid_out, true);
460
461
462
463
464
465
466
		grib_handle_delete(h);
	} else {
		cleanup();
		throw IOException("No grid found in file \""+filename+"\"", AT);
	}

	cleanup();
467
468
}

469
void GRIBIO::indexFile(const std::string& filename)
470
{
471
472
473
474
475
476
477
	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);
	}

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

	//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);
	}

508
}
509

510
void GRIBIO::read2DGrid(Grid2DObject& grid_out, const MeteoGrids::Parameters& parameter, const Date& date)
511
{
512
513
	Date UTC_date = date;
	UTC_date.setTimeZone(tz_in);
514

515
	const std::string filename = grid2dpath_in+"/"+grid2d_prefix+UTC_date.toString(Date::NUM).substr(0,10)+grid2d_ext;
516
517

	read2DGrid(filename, grid_out, parameter, UTC_date);
518
519
}

520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
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++) {
536
				VW(ii,jj) = sqrt( Optim::pow2(U(ii,jj)) + Optim::pow2(V(ii,jj)) );
537
				DW(ii,jj) = fmod( atan2( U(ii,jj), V(ii,jj) ) * Cst::to_deg + 360. + bearing_offset, 360.); // turn into degrees [0;360)
538
539
540
541
542
543
544
			}
		}
	}

	wind_date = date;
}

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

552
	//Basic meteo parameters
553
	if(parameter==MeteoGrids::P) read2DGrid_indexed(1.2, 1, 0, date, grid_out); //PS
554
	if(parameter==MeteoGrids::TA) read2DGrid_indexed(11.2, 105, 2, date, grid_out); //T_2M
555
	if(parameter==MeteoGrids::RH) {
556
		if(!read2DGrid_indexed(52.2, 105, 2, date, grid_out)) { //RELHUM_2M
557
			Grid2DObject ta;
558
559
			read2DGrid_indexed(11.2, 105, 2, date, ta); //T_2M
			read2DGrid_indexed(17.2, 105, 2, date, grid_out); //TD_2M
560
561
562
563
564
565
566
			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);
				}
			}
		}
	}
567
568
	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
569
570

	//hydrological parameters
571
572
573
574
	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) {
575
576
577
578
579
580
		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;
		}
581
	}
582
583

	//radiation parameters
584
585
586
587
	if(parameter==MeteoGrids::ALB) {
		read2DGrid_indexed(84.2, 1, 0, date, grid_out); //ALB_RAD
		grid_out.grid2D /= 100.;
	}
588
589
590
591
592
	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
	}
593
	/*if(parameter==MeteoGrids::CLD) { //cloudiness
594
595
596
597
		if(read2DGrid_indexed(74.2, 1, 0, date, grid_out)) //CLCM
		grid_out.grid2D /= 100.;
	}*/

598
	if(parameter==MeteoGrids::ISWR) {
599
600
601
602
		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
603
			Grid2DObject diff;
604
605
			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
606
607
608
609
610
			grid_out.grid2D += diff.grid2D;
		}
	}

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

625
	//Wind parameters
626
	if(parameter==MeteoGrids::VW_MAX) read2DGrid_indexed(187.201, 105, 10, date, grid_out); //VMAX_10M 10m
627
628
629
630
631
632
	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++) {
633
				grid_out(ii,jj) = VW(ii,jj)*sin(DW(ii,jj)*Cst::to_rad);
634
635
636
			}
		}
	}
637
638
639
640
	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++) {
641
				grid_out(ii,jj) = VW(ii,jj)*cos(DW(ii,jj)*Cst::to_rad);
642
643
644
			}
		}
	}
645
646
647
648
649
650
651
652
	if(parameter==MeteoGrids::DW) {
		readWind(filename, date);
		grid_out = DW;
	}
	if(parameter==MeteoGrids::VW) {
		readWind(filename, date);
		grid_out = VW;
	}
653
654
655
656
657
658
659
660
661

	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
662
663
	/*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
664
		Grid2DObject Z0;
665
		if(read2DGrid_indexed(83.2, 1, 0, date, Z0)) { //Z0
666
667
			for(unsigned int jj=0; jj<grid_out.nrows; jj++) {
				for(unsigned int ii=0; ii<grid_out.ncols; ii++) {
668
					grid_out(ii,jj) = Atmosphere::windLogProfile(grid_out(ii,jj), 10., 7.5, Z0(ii,jj));
669
670
				}
			}
671
		} else {
672
			const double wind_factor = Atmosphere::windLogProfile(1., 10., 7.5, 0.03);
673
			grid_out.grid2D *= wind_factor;
674
		}
675
	}*/
676
677
}

678
void GRIBIO::readDEM(DEMObject& dem_out)
679
{
680
	const Date d; //ie: undef. This will be caught when reading the GRIB file
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
	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();
	}
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
}

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);
}

729
730
731
void GRIBIO::scanMeteoPath()
{
	std::list<std::string> dirlist;
732
	IOUtils::readDirectory(meteopath_in, dirlist, meteo_ext);
733
734
735
736
737
738
739
740
	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;
741
		IOUtils::convertString(date, filename.substr(spos,10), tz_in);
742
743
744
745
746
747
748
		std::pair<Date,std::string> tmp(date, filename);

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

749
750
void GRIBIO::readMeteoData(const Date& dateStart, const Date& dateEnd,
                             std::vector< std::vector<MeteoData> >& vecMeteo,
751
752
                             const size_t&)
{
753
754
755
756
757
	if(!meteo_initialized) {
		readStations();
		scanMeteoPath();
		meteo_initialized=true;
	}
758

759
	vecMeteo.clear();
760

761
762
	double *lats = (double*)malloc(vecPts.size()*sizeof(double));
	double *lons = (double*)malloc(vecPts.size()*sizeof(double));
763
764
765
766
767
768
769
770
771
772
773
774
775
776
	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)
777
778

	try {
779
780
781
782
		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;
783
784
785

			if(!indexed || idx_filename!=filename) {
				cleanup();
786
				indexFile(filename); //this will also read geolocalization
787
788
			}
			if(!meta_ok) {
789
				if(readMeteoMeta(vecPts, meta, lats, lons)==false) {
790
791
792
793
					//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));
794
					readMeteoMeta(vecPts, meta, lats, lons);
795
				}
796
				vecMeteo.insert(vecMeteo.begin(), vecPts.size(), std::vector<MeteoData>()); //allocation for the vectors now that we know how many true stations we have
797
				meta_ok=true;
798
799
800
801
802
803
804
805
806
			}

			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);
807
		cleanup();
808
		throw;
809
810
	}

811
	free(lats); free(lons);
812
813
}

814
815
816
817
818
819
820
821
822
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++) {
823
824
825
			if(lat==lats[jj] && lon==lons[jj]) {
				deletions.push_back(jj);
			}
826
827
828
829
830
831
		}
	}

	//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];
832
		vecPts.erase(vecPts.begin()+index);
833
834
835
836
837
838
	}

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

839
bool GRIBIO::readMeteoMeta(std::vector<Coords>& vecPts, std::vector<StationData> &stations, double *lats, double *lons)
840
841
842
{//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();

843
844
845
846
	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;
847
848
849
850
851
	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);
	}
852

853
	const long npoints = vecPts.size();
854
	double latitudeOfSouthernPole, longitudeOfSouthernPole;
855
856
	GRIB_CHECK(grib_get_double(h,"latitudeOfSouthernPoleInDegrees",&latitudeOfSouthernPole),0);
	GRIB_CHECK(grib_get_double(h,"longitudeOfSouthernPoleInDegrees",&longitudeOfSouthernPole),0);
857
858
	latitudeOfNorthernPole = -latitudeOfSouthernPole;
	longitudeOfNorthernPole = longitudeOfSouthernPole+180.;
859

860
861
	long Ni;
	GRIB_CHECK(grib_get_long(h,"Ni",&Ni),0);
862

863
864
	//build GRIB local coordinates for the points
	for(unsigned int ii=0; ii<(unsigned)npoints; ii++) {
865
		Coords::trueLatLonToRotated(latitudeOfNorthernPole, longitudeOfNorthernPole, vecPts[ii].getLat(), vecPts[ii].getLon(), lats[ii], lons[ii]);
866
	}
867

868
869
870
871
872
873
874
875
876
877
878
	//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);
	}
879

880
881
	//remove potential duplicates
	if(removeDuplicatePoints(vecPts, outlats, outlons)==true) {
882
883
		free(outlats); free(outlons); free(values); free(distances); free(indexes);
		grib_handle_delete(h);
884
		return false;
885
	}
886
887
888
889
890

	//fill metadata
	for(unsigned int ii=0; ii<(unsigned)npoints; ii++) {
		StationData sd;
		sd.position.setProj(coordin, coordinparam);
891
		double true_lat, true_lon;
892
		Coords::rotatedToTrueLatLon(latitudeOfNorthernPole, longitudeOfNorthernPole, outlats[ii], outlons[ii], true_lat, true_lon);
893
		sd.position.setLatLon(true_lat, true_lon, values[ii]);
894
895
896
897
898
899
900
901
902
903
904
905
		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;
906
907
908
909
910
911
}

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);
912
913
914
915
916
917
918
919
920

	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);
		}

921
922
923
924
925
926
927
928
		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);

929
930
		long level=0;
		if(i_level!=0) GRIB_CHECK(grib_get_long(h,"level", &level),0);
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
		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);
948
				grib_handle_delete(h);
949
				return true;
950
951
952
953
954
			}
		}
		grib_handle_delete(h);
	}
	return false;
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
987
988
989
990
991
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
992
993
	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
994
995
996
997
998
999
1000
	   && 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
For faster browsing, not all history is shown. View entire blame