WSL/SLF GitLab Repository

GRIBIO.cc 41.8 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
#include <meteoio/ResamplingAlgorithms2D.h>
21
22
23
#include <meteoio/meteoLaws/Atmosphere.h>
#include <meteoio/meteoLaws/Meteoconst.h> //for PI
#include <meteoio/dataClasses/DEMObject.h>
24
#include <meteoio/MathOptim.h>
25

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

using namespace std;

namespace mio {
/**
35
 * @page gribio GRIBIO
36
 * @section gribio_format Format and limitations
37
 * This plugin reads GRIB (https://en.wikipedia.org/wiki/GRIB) files as produced by meteorological models.
38
 * Being based on GRIB API (https://software.ecmwf.int/wiki/display/GRIB/Home), 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)).
39
 * 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://rda.ucar.edu/docs/formats/grib/gribdoc/params.html) and levels
40
 * (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 .
41
 *
42
 * Several assumptions/approximations are held/made when reading grids:
43
 * - 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".
44
45
46
47
 * - 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.
 *
48
 * 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 .
49
50
 *
 * As a side note, when calling read2DGrid(grid, filename), it will returns the first grid that is found.
51
 *
52
 * @section gribio_cosmo_partners COSMO Group
53
54
55
56
57
58
59
60
61
62
63
64
65
66
 * 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
 *
67
68
 * @section gribio_units Units
 * As specified by WMO.
69
 *
70
 * @section gribio_files_naming Files naming scheme
Mathias Bavay's avatar
Mathias Bavay committed
71
72
73
 * When a grid is read by providing the filename to open, any file name will obviously work. Otherwise, the file names have to follow the pattern:\n
 * {GRID2DPREFIX}{ISO8601 numerical UTC date}{GRID2DEXT}\n
 * By default, GRID2DPREFIX is empty and GRID2DEXT is ".grb". This means that by default, a grib file containing data for 2013-10-15T12:00 would be:
74
75
 * "201310151200.grb". Since the grid extension contains the ".", it is possible to use it for a file name suffix as well. For example, to read files
 * named like "lfff201310151200_z" one would provide GRID2DPREFIX="lfff" and GRID2DEXT="_z".
Mathias Bavay's avatar
Mathias Bavay committed
76
 *
77
 * @section gribio_keywords Keywords
78
 * This plugin uses the following keywords:
79
80
81
 * - COORDSYS: coordinate system (see Coords)
 * - COORDPARAM: extra coordinates parameters (see Coords)
 * - GRID2DPATH: path where to find the grids
82
 * - GRID2DPREFIX: prefix to append when generating a file name for reading (ie: something like "laf" for Cosmo-Analysis-full domain), optional
Mathias Bavay's avatar
Mathias Bavay committed
83
 * - GRID2DEXT: grib file extension, or <i>none</i> for no file extension (default: ".grb")
84
85
 * - 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)
86
 * - METEOPATH: path where to find the grids for extracting time series at the given points
Mathias Bavay's avatar
Mathias Bavay committed
87
 * - METEOEXT: file extension, or <i>none</i> for no file extension (default: ".grb")
88
 * - STATION#: coordinates for virtual stations (if using GRIB as METEO plugin). Each station is given by its coordinates and the closest
89
90
 * 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.
91
 *
92
93
94
 */

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

98
99
100
101
102
GRIBIO::GRIBIO(const std::string& configfile)
        : cfg(configfile), grid2dpath_in(), meteopath_in(), vecPts(), cache_meteo_files(),
          meteo_ext(default_ext), grid2d_ext(default_ext), grid2d_prefix(), idx_filename(), coordin(), coordinparam(),
          VW(), DW(), wind_date(), llcorner(), fp(NULL), idx(NULL),
          latitudeOfNorthernPole(IOUtils::nodata), longitudeOfNorthernPole(IOUtils::nodata), bearing_offset(IOUtils::nodata),
103
          cellsize_x(IOUtils::nodata), cellsize_y(IOUtils::nodata), factor_x(IOUtils::nodata), factor_y(IOUtils::nodata),
104
          indexed(false), meteo_initialized(false), update_dem(false)
105
106
107
108
{
	setOptions();
}

109
110
111
112
113
GRIBIO::GRIBIO(const Config& cfgreader)
        : cfg(cfgreader), grid2dpath_in(), meteopath_in(), vecPts(), cache_meteo_files(),
          meteo_ext(default_ext), grid2d_ext(default_ext), grid2d_prefix(), idx_filename(), coordin(), coordinparam(),
          VW(), DW(), wind_date(), llcorner(), fp(NULL), idx(NULL),
          latitudeOfNorthernPole(IOUtils::nodata), longitudeOfNorthernPole(IOUtils::nodata), bearing_offset(IOUtils::nodata),
114
          cellsize_x(IOUtils::nodata), cellsize_y(IOUtils::nodata), factor_x(IOUtils::nodata), factor_y(IOUtils::nodata),
115
          indexed(false), meteo_initialized(false), update_dem(false)
116
117
{
	setOptions();
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
}

GRIBIO& GRIBIO::operator=(const GRIBIO& source) {
	if(this != &source) {
		fp = NULL;
		idx = NULL;
		grid2dpath_in = source.grid2dpath_in;
		meteopath_in = source.meteopath_in;
		vecPts = source.vecPts;
		cache_meteo_files = source.cache_meteo_files;
		meteo_ext = source.meteo_ext;
		grid2d_ext = source.grid2d_ext;
		grid2d_prefix = source.grid2d_prefix;
		idx_filename = source.idx_filename;
		coordin = source.coordin;
		coordinparam = source.coordinparam;
		VW = source.VW;
		DW = source.DW;
		wind_date = source.wind_date;
		llcorner = source.llcorner;
		latitudeOfNorthernPole = source.latitudeOfNorthernPole;
		longitudeOfNorthernPole = source.longitudeOfNorthernPole;
		bearing_offset = source.bearing_offset;
		cellsize_x = source.cellsize_x;
		cellsize_y = source.cellsize_y;
143
144
		factor_x = source.factor_x;
		factor_y = source.factor_y;
145
146
147
148
149
		indexed = source.indexed;
		meteo_initialized = source.meteo_initialized;
		update_dem = source.update_dem;
	}
	return *this;
150
151
152
153
}

GRIBIO::~GRIBIO() throw()
{
154
	cleanup();
155
156
157
158
}

void GRIBIO::setOptions()
{
159
	std::string coordout, coordoutparam;
160
	IOUtils::getProjectionParameters(cfg, coordin, coordinparam, coordout, coordoutparam);
161

162
	const string tmp = cfg.get("GRID2D", "Input", IOUtils::nothrow);
163
164
	if (tmp == "GRIB") { //keep it synchronized with IOHandler.cc for plugin mapping!!
		cfg.getValue("GRID2DPATH", "Input", grid2dpath_in);
165
		cfg.getValue("GRIB_DEM_UPDATE", "Input", update_dem, IOUtils::nothrow);
166
	}
167
	cfg.getValue("GRID2DPREFIX", "Input", grid2d_prefix, IOUtils::nothrow);
168

169
	cfg.getValue("METEOEXT", "Input", meteo_ext, IOUtils::nothrow);
170
	if(meteo_ext=="none") meteo_ext.clear();
171

172
	cfg.getValue("GRID2DEXT", "Input", grid2d_ext, IOUtils::nothrow);
173
	if(grid2d_ext=="none") grid2d_ext.clear();
174
}
175

176
void GRIBIO::readStations(std::vector<Coords> &vecPoints)
177
178
{
	cfg.getValue("METEOPATH", "Input", meteopath_in);
179

180
181
182
183
184
185
	std::vector<std::string> vecStation;
	cfg.getValues("STATION", "INPUT", vecStation);
	for(size_t ii=0; ii<vecStation.size(); ii++) {
		Coords tmp(coordin, coordinparam, vecStation[ii]);
		if(!tmp.isNodata())
			vecPoints.push_back( tmp );
186

187
188
		std::cerr <<  "\tRead virtual station " << vecPoints.back().printLatLon() << "\n";
	}
189
190
}

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

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

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

	grib_keys_iterator_delete(kiter);
}

216
void  GRIBIO::listFields(const std::string& filename)
217
{
218
	grib_handle* h=NULL;
219
	int err=0;
220

221
	//loop over the messages
222
223
224
225
226
227
	while((h = grib_handle_new_from_file(0,fp,&err)) != NULL) {
		if(!h) {
			cleanup();
			throw IOException("Unable to create grib handle for \""+filename+"\"", AT);
		}

228
229
230
231
232
233
234
235
		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);
236
		long levelType;
237
238
239
		GRIB_CHECK(grib_get_long(h,"indicatorOfTypeOfLevel", &levelType),0);
		long level=0;
		if(levelType!=1) GRIB_CHECK(grib_get_long(h,"level", &level),0);
240
		std::cerr << marsParam << " " << shortname << " " << name << " type " << levelType << " level " << level << "\n";
241
		grib_handle_delete(h);
242
243
	}
}
244

245
void GRIBIO::getDate(grib_handle* h, Date &base, double &d1, double &d2) {
246
247
248
249
	long dataDate, dataTime;
	GRIB_CHECK(grib_get_long(h,"dataDate",&dataDate),0);
	GRIB_CHECK(grib_get_long(h,"dataTime",&dataTime),0);

250
251
	const int year=static_cast<int>(dataDate/10000), month=static_cast<int>(dataDate/100-year*100), day=static_cast<int>(dataDate-month*100-year*10000);
	const int hour=static_cast<int>(dataTime/100), minutes=static_cast<int>(dataTime-hour*100);
252
253
254
255
256
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
	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:
290
			std::ostringstream ss;
291
292
293
294
			ss << "GRIB file using stepUnits=" << stepUnits << ", which is not supported";
			throw InvalidFormatException(ss.str(), AT);
	}

295
296
	d1 = static_cast<double>(startStep*step_units);
	d2 = static_cast<double>(endStep*step_units);
297
298
}

299
Coords GRIBIO::getGeolocalization(grib_handle* h, double &cell_x, double &cell_y)
300
{
301
	//getting transformation parameters
302
303
	double angleOfRotationInDegrees = 0.; //if the key was not found, we consider there is no rotation
	grib_get_double(h,"angleOfRotationInDegrees",&angleOfRotationInDegrees);
304
305
	if(angleOfRotationInDegrees!=0.) {
		throw InvalidArgumentException("Rotated grids not supported!", AT);
306
	}
307
308
309
	double latitudeOfSouthernPole = -90, longitudeOfSouthernPole = -180;
	grib_get_double(h,"latitudeOfSouthernPoleInDegrees",&latitudeOfSouthernPole);
	grib_get_double(h,"longitudeOfSouthernPoleInDegrees",&longitudeOfSouthernPole);
310
311
	latitudeOfNorthernPole = -latitudeOfSouthernPole;
	longitudeOfNorthernPole = longitudeOfSouthernPole+180.;
312

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

	//determining cell size
326
327
328
	long Ni, Nj;
	GRIB_CHECK(grib_get_long(h,"Ni",&Ni),0);
	GRIB_CHECK(grib_get_long(h,"Nj",&Nj),0);
329
	double bearing;
330
331
332
333
334
335
336
337
	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)
338

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

	return cntr; //this is now the ll corner
346
347
}

348
void GRIBIO::read2Dlevel(grib_handle* h, Grid2DObject& grid_out, const bool& read_geolocalization)
349
{
350
351
352
353
354
355
356
	long Ni, Nj;
	GRIB_CHECK(grib_get_long(h,"Ni",&Ni),0);
	GRIB_CHECK(grib_get_long(h,"Nj",&Nj),0);

	size_t values_len= 0;
	GRIB_CHECK(grib_get_size(h,"values",&values_len),0);
	if(values_len!=(unsigned)(Ni*Nj)) {
357
		ostringstream ss;
358
359
360
361
		ss << "Declaring grid of size " << Ni << "x" << Nj << "=" << Ni*Nj << " ";
		ss << "but containing " << values_len << " values. This is inconsistent!";
		throw InvalidArgumentException(ss.str(), AT);
	}
362
	double *values = (double*)calloc(values_len, sizeof(double));
363
364

	GRIB_CHECK(grib_get_double_array(h,"values",values,&values_len),0);
365

366
	if(read_geolocalization) llcorner = getGeolocalization(h, cellsize_x, cellsize_y);
367
368
369
370
	grid_out.set(static_cast<size_t>(Ni), static_cast<size_t>(Nj), .5*(cellsize_x+cellsize_y), llcorner);
	size_t i=0;
	for(size_t jj=0; jj<(unsigned)Nj; jj++) {
		for(size_t ii=0; ii<(unsigned)Ni; ii++)
371
			grid_out(ii,jj) = values[i++];
372
	}
373
374
375
376
	//cells were not square, we have to resample
	if (factor_x!=IOUtils::nodata && factor_y!=IOUtils::nodata) {
		grid_out.grid2D = ResamplingAlgorithms2D::BilinearResampling(grid_out.grid2D, factor_x, factor_y);
	}
377
378
	free(values);
}
379

380
bool GRIBIO::read2DGrid_indexed(const double& in_marsParam, const long& i_levelType, const long& i_level, const Date i_date, Grid2DObject& grid_out)
381
{
382
383
384
	GRIB_CHECK(grib_index_select_double(idx,"marsParam",in_marsParam),0);
	GRIB_CHECK(grib_index_select_long(idx,"indicatorOfTypeOfLevel", i_levelType),0);

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

393
394
395
396
397
398
399
400
		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);

401
		long level=0;
402
		if(i_level!=0) GRIB_CHECK(grib_get_long(h,"level", &level),0);
403
404
405
406
407
408
409
410
411
412
		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;
			}
413
		}
414
		grib_handle_delete(h);
415
	}
416
	return false;
417
418
}

419
void GRIBIO::read2DGrid(Grid2DObject& grid_out, const std::string& i_name)
420
{
421
422
423
	const std::string filename = grid2dpath_in+"/"+i_name;
	fp = fopen(filename.c_str(),"r");
	if(fp==NULL) {
424
		ostringstream ss;
425
		ss << "Error opening file \"" << filename << "\", possible reason: " << strerror(errno);
426
427
428
429
430
431
432
433
434
435
436
		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);
		}

437
		read2Dlevel(h, grid_out, true);
438
439
440
441
442
443
444
		grib_handle_delete(h);
	} else {
		cleanup();
		throw IOException("No grid found in file \""+filename+"\"", AT);
	}

	cleanup();
445
446
}

447
void GRIBIO::indexFile(const std::string& filename)
448
{
449
450
	fp = fopen(filename.c_str(),"r");
	if(fp==NULL) {
451
		ostringstream ss;
452
		ss << "Error opening file \"" << filename << "\", possible reason: " << strerror(errno);
453
454
455
		throw FileAccessException(ss.str(), AT);
	}

456
	int err=0;
457
	const std::string keys("marsParam:d,indicatorOfTypeOfLevel:l"); //indexing keys
458
459
460
461
	char *c_filename = (char *)filename.c_str();
	idx = grib_index_new_from_file(0, c_filename, keys.c_str(), &err);
	if(err!=0) {
		cleanup();
462
		throw IOException("Failed to index GRIB file \""+filename+"\". Is it a valid GRIB file?", AT);
463
464
465
	}
	indexed=true;
	idx_filename=filename;
466
467
468
469
470
471
472
473
474
475

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

476
		llcorner = getGeolocalization(h, cellsize_x, cellsize_y); //this sets llcorner, cellsize and bearing_offset
477
		if( fabs(cellsize_x-cellsize_y)/cellsize_x > 1./100.) {
478
479
480
481
			// round to 1cm precision for numerical stability
			const double cellsize =  Optim::round( std::min(cellsize_x, cellsize_y) * 100. ) / 100.;
			factor_x = cellsize_x / cellsize;
			factor_y = cellsize_y / cellsize;
482
483
484
485
486
487
488
		}
		grib_handle_delete(h);
	} else {
		cleanup();
		throw IOException("No grid found in file \""+filename+"\"", AT);
	}

489
}
490

491
void GRIBIO::read2DGrid(Grid2DObject& grid_out, const MeteoGrids::Parameters& parameter, const Date& date)
492
{
493
494
	Date UTC_date = date;
	UTC_date.setTimeZone(tz_in);
495

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

	read2DGrid(filename, grid_out, parameter, UTC_date);
499
500
}

501
502
503
504
505
506
507
508
509
510
511
512
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

513
514
515
516
		VW.set(U.getNx(), U.getNy(), U.cellsize, U.llcorner);
		DW.set(U.getNx(), U.getNy(), U.cellsize, U.llcorner);
		for(size_t jj=0; jj<VW.getNy(); jj++) {
			for(size_t ii=0; ii<VW.getNx(); ii++) {
517
				VW(ii,jj) = sqrt( Optim::pow2(U(ii,jj)) + Optim::pow2(V(ii,jj)) );
518
				DW(ii,jj) = fmod( atan2( U(ii,jj), V(ii,jj) ) * Cst::to_deg + 360. + bearing_offset, 360.); // turn into degrees [0;360)
519
520
521
522
523
524
525
			}
		}
	}

	wind_date = date;
}

526
void GRIBIO::read2DGrid(const std::string& filename, Grid2DObject& grid_out, const MeteoGrids::Parameters& parameter, const Date& date)
527
{ //Parameters should be read in table 2 if available since this table is the one standardized by WMO
528
	if(!indexed || idx_filename!=filename) {
529
530
531
532
		cleanup();
		indexFile(filename);
	}

533
	//Basic meteo parameters
534
	if(parameter==MeteoGrids::P) read2DGrid_indexed(1.2, 1, 0, date, grid_out); //PS
535
	if(parameter==MeteoGrids::TA) read2DGrid_indexed(11.2, 105, 2, date, grid_out); //T_2M
536
	if(parameter==MeteoGrids::RH) {
537
		if(!read2DGrid_indexed(52.2, 105, 2, date, grid_out)) { //RELHUM_2M
538
			Grid2DObject ta;
539
540
			read2DGrid_indexed(11.2, 105, 2, date, ta); //T_2M
			read2DGrid_indexed(17.2, 105, 2, date, grid_out); //TD_2M
541
542
			for(size_t jj=0; jj<grid_out.getNy(); jj++) {
				for(size_t ii=0; ii<grid_out.getNx(); ii++) {
543
544
545
546
547
					grid_out(ii,jj) = Atmosphere::DewPointtoRh(grid_out(ii,jj), ta(ii,jj), true);
				}
			}
		}
	}
548
549
	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
550
551

	//hydrological parameters
552
553
554
555
	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) {
556
557
558
559
560
561
		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;
		}
562
	}
563
564

	//radiation parameters
565
566
567
568
	if(parameter==MeteoGrids::ALB) {
		read2DGrid_indexed(84.2, 1, 0, date, grid_out); //ALB_RAD
		grid_out.grid2D /= 100.;
	}
569
570
571
572
573
	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
	}
574
	/*if(parameter==MeteoGrids::CLD) { //cloudiness
575
576
577
578
		if(read2DGrid_indexed(74.2, 1, 0, date, grid_out)) //CLCM
		grid_out.grid2D /= 100.;
	}*/

579
	if(parameter==MeteoGrids::ISWR) {
580
581
		if(read2DGrid_indexed(116.2, 1, 0, date, grid_out)) { //short wave
			grid_out.grid2D *= -1.;
582
		} else if(!read2DGrid_indexed(111.250, 1, 0, date, grid_out)) { //GLOB
583
			Grid2DObject diff;
584
585
			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
586
587
588
589
590
			grid_out.grid2D += diff.grid2D;
		}
	}

	//DEM parameters
591
	if(parameter==MeteoGrids::DEM) read2DGrid_indexed(8.2, 1, 0, date, grid_out); //HSURF
592
593
	if(parameter==MeteoGrids::SLOPE) {
		read2DGrid_indexed(98.202, 1, 0, date, grid_out); //SLO_ANG
594
		grid_out.grid2D *= Cst::to_deg;
595
596
597
	}
	if(parameter==MeteoGrids::AZI) {
		read2DGrid_indexed(99.202, 1, 0, date, grid_out); //SLO_ASP
598
599
		for(size_t jj=0; jj<grid_out.getNy(); jj++) {
			for(size_t ii=0; ii<grid_out.getNx(); ii++) {
600
				grid_out(ii,jj) = fmod( grid_out(ii,jj)*Cst::to_deg + 360. + bearing_offset, 360.); // turn into degrees [0;360)
601
602
			}
		}
603
	}
604

605
	//Wind parameters
606
	if(parameter==MeteoGrids::VW_MAX) read2DGrid_indexed(187.201, 105, 10, date, grid_out); //VMAX_10M 10m
607
608
609
610
	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);
611
612
		for(size_t jj=0; jj<grid_out.getNy(); jj++) {
			for(size_t ii=0; ii<grid_out.getNx(); ii++) {
613
				grid_out(ii,jj) = VW(ii,jj)*sin(DW(ii,jj)*Cst::to_rad);
614
615
616
			}
		}
	}
617
618
	if(parameter==MeteoGrids::V) {
		readWind(filename, date);
619
620
		for(size_t jj=0; jj<grid_out.getNy(); jj++) {
			for(size_t ii=0; ii<grid_out.getNx(); ii++) {
621
				grid_out(ii,jj) = VW(ii,jj)*cos(DW(ii,jj)*Cst::to_rad);
622
623
624
			}
		}
	}
625
626
627
628
629
630
631
632
	if(parameter==MeteoGrids::DW) {
		readWind(filename, date);
		grid_out = DW;
	}
	if(parameter==MeteoGrids::VW) {
		readWind(filename, date);
		grid_out = VW;
	}
633
634

	if(grid_out.isEmpty()) {
635
		ostringstream ss;
636
637
638
639
640
641
		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
642
643
	/*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
644
		Grid2DObject Z0;
645
		if(read2DGrid_indexed(83.2, 1, 0, date, Z0)) { //Z0
Mathias Bavay's avatar
Mathias Bavay committed
646
647
			for(size_t jj=0; jj<grid_out.getNy(); jj++) {
				for(size_t ii=0; ii<grid_out.getNx(); ii++) {
648
					grid_out(ii,jj) = Atmosphere::windLogProfile(grid_out(ii,jj), 10., 7.5, Z0(ii,jj));
649
650
				}
			}
651
		} else {
652
			const double wind_factor = Atmosphere::windLogProfile(1., 10., 7.5, 0.03);
653
			grid_out.grid2D *= wind_factor;
654
		}
655
	}*/
656
657
}

658
void GRIBIO::readDEM(DEMObject& dem_out)
659
{
660
	const Date d(2012,9,19,0,0,0); //ie: undef. This will be caught when reading the GRIB file
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
	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();
	}
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
}

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

709
710
711
void GRIBIO::scanMeteoPath()
{
	std::list<std::string> dirlist;
712
	IOUtils::readDirectory(meteopath_in, dirlist, meteo_ext);
713
714
715
	dirlist.sort();

	//Check date in every filename and cache it
716
	std::list<std::string>::const_iterator it = dirlist.begin();
717
718
	while ((it != dirlist.end())) {
		const std::string& filename = *it;
719
		const std::string::size_type spos = filename.find_first_of("0123456789");
720
		Date date;
721
		IOUtils::convertString(date, filename.substr(spos,10), tz_in);
722
		const std::pair<Date,std::string> tmp(date, filename);
723
724
725
726
727
728

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

729
730
void GRIBIO::readMeteoData(const Date& dateStart, const Date& dateEnd,
                             std::vector< std::vector<MeteoData> >& vecMeteo,
731
732
                             const size_t&)
{
733
	if(!meteo_initialized) {
734
		readStations(vecPts);
735
736
737
		scanMeteoPath();
		meteo_initialized=true;
	}
738

739
	vecMeteo.clear();
740

741
742
	double *lats = (double*)malloc(vecPts.size()*sizeof(double));
	double *lons = (double*)malloc(vecPts.size()*sizeof(double));
743
744
745
746
	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
747
	size_t idx_start;
748
	bool start_found=false;
Thomas Egger's avatar
Thomas Egger committed
749
	for (idx_start=0; idx_start<cache_meteo_files.size(); idx_start++) {
750
751
752
753
754
		if(dateStart<cache_meteo_files[idx_start].first) {
			start_found=true;
			break;
		}
	}
Thomas Egger's avatar
Thomas Egger committed
755
756
757
758
759
760
761

	if (start_found==false) {
		free(lats); free(lons);
		return;
	}

	if (idx_start>0) idx_start--; //start with first element before dateStart (useful for resampling)
762
763

	try {
764
		for(size_t ii=idx_start; ii<cache_meteo_files.size(); ii++) {
765
766
767
			const Date& date = cache_meteo_files[ii].first;
			if(date>dateEnd) break;
			const std::string filename = meteopath_in+"/"+cache_meteo_files[ii].second;
768
769
770

			if(!indexed || idx_filename!=filename) {
				cleanup();
771
				indexFile(filename); //this will also read geolocalization
772
773
			}
			if(!meta_ok) {
774
				if(readMeteoMeta(vecPts, meta, lats, lons)==false) {
775
776
777
778
					//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));
779
					readMeteoMeta(vecPts, meta, lats, lons);
780
				}
781
				vecMeteo.insert(vecMeteo.begin(), vecPts.size(), std::vector<MeteoData>()); //allocation for the vectors now that we know how many true stations we have
782
				meta_ok=true;
783
784
785
786
			}

			std::vector<MeteoData> Meteo;
			readMeteoStep(meta, lats, lons, date, Meteo);
787
			for(size_t jj=0; jj<vecPts.size(); jj++)
788
789
790
791
				vecMeteo[jj].push_back(Meteo[jj]);
		}
	} catch(...) {
		free(lats); free(lons);
792
		cleanup();
793
		throw;
794
795
	}

796
	free(lats); free(lons);
797
798
}

799
bool GRIBIO::removeDuplicatePoints(std::vector<Coords>& vecPoints, double *lats, double *lons)
800
{ //remove potential duplicates. Returns true if some have been removed
801
	const size_t npoints = vecPoints.size();
802
803
	std::vector<size_t> deletions;
	deletions.reserve(npoints);
804
	for(size_t ii=0; ii<npoints; ii++) {
805
806
		const double lat = lats[ii];
		const double lon = lons[ii];
807
		for(size_t jj=ii+1; jj<npoints; jj++) {
808
809
810
			if(lat==lats[jj] && lon==lons[jj]) {
				deletions.push_back(jj);
			}
811
812
813
814
		}
	}

	//we need to erase from the end in order to keep the index unchanged...
815
816
	for(size_t ii=deletions.size(); ii>0; ii--) {
		const size_t index=deletions[ii-1];
817
		vecPoints.erase(vecPoints.begin()+index);
818
819
	}

820
	if(!deletions.empty()) return true;
821
822
823
	return false;
}

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

828
829
830
831
	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;
832
833
834
835
836
	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);
	}
837

838
	const size_t npoints = vecPoints.size();
839
	double latitudeOfSouthernPole, longitudeOfSouthernPole;
840
841
	GRIB_CHECK(grib_get_double(h,"latitudeOfSouthernPoleInDegrees",&latitudeOfSouthernPole),0);
	GRIB_CHECK(grib_get_double(h,"longitudeOfSouthernPoleInDegrees",&longitudeOfSouthernPole),0);
842
843
	latitudeOfNorthernPole = -latitudeOfSouthernPole;
	longitudeOfNorthernPole = longitudeOfSouthernPole+180.;
844

845
846
	long Ni;
	GRIB_CHECK(grib_get_long(h,"Ni",&Ni),0);
847

848
	//build GRIB local coordinates for the points
849
	for(size_t ii=0; ii<npoints; ii++) {
850
		Coords::trueLatLonToRotated(latitudeOfNorthernPole, longitudeOfNorthernPole, vecPoints[ii].getLat(), vecPoints[ii].getLon(), lats[ii], lons[ii]);
851
	}
852

853
854
855
856
857
858
859
860
861
862
863
	//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);
	}
864

865
	//remove potential duplicates
866
	if(removeDuplicatePoints(vecPoints, outlats, outlons)==true) {
867
868
		free(outlats); free(outlons); free(values); free(distances); free(indexes);
		grib_handle_delete(h);
869
		return false;
870
	}
871
872

	//fill metadata
873
	for(size_t ii=0; ii<npoints; ii++) {
874
875
		StationData sd;
		sd.position.setProj(coordin, coordinparam);
876
		double true_lat, true_lon;
877
		Coords::rotatedToTrueLatLon(latitudeOfNorthernPole, longitudeOfNorthernPole, outlats[ii], outlons[ii], true_lat, true_lon);
878
		sd.position.setLatLon(true_lat, true_lon, values[ii]);
879
		ostringstream ss;
880
881
		ss << "Point_" << indexes[ii];
		sd.stationID=ss.str();
882
		ostringstream ss2;
883
884
885
886
887
888
889
890
		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;
891
892
}

893
bool GRIBIO::readMeteoValues(const double& marsParam, const long& levelType, const long& i_level, const Date& i_date, const size_t& npoints, double *lats, double *lons, double *values)
894
895
896
{
	GRIB_CHECK(grib_index_select_double(idx,"marsParam",marsParam),0);
	GRIB_CHECK(grib_index_select_long(idx,"indicatorOfTypeOfLevel", levelType),0);
897
898
899
900
901
902
903
904
905

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

906
907
908
909
910
911
912
913
		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);

914
915
		long level=0;
		if(i_level!=0) GRIB_CHECK(grib_get_long(h,"level", &level),0);
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
		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);
933
				grib_handle_delete(h);
934
				return true;
935
936
937
938
939
			}
		}
		grib_handle_delete(h);
	}
	return false;
940
941
}

942
943
void GRIBIO::fillMeteo(double *values, const MeteoData::Parameters& param, const size_t& npoints, std::vector<MeteoData> &Meteo) {
	for(size_t ii=0; ii<npoints; ii++) {
944
945
946
947
948
949
		Meteo[ii](param) = values[ii];
	}
}

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