WSL/SLF GitLab Repository

GRIBIO.cc 42.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
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 (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)).
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://rda.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
 * Several assumptions/approximations are held/made when reading grids:
41
 * - 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".
42
43
44
45
 * - 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.
 *
46
 * 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 .
47
48
 *
 * 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
 *
Mathias Bavay's avatar
Mathias Bavay committed
68
69
70
71
 * @section files_naming Files naming scheme
 * 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:
72
73
 * "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
74
 *
75
 * @section gribio_keywords Keywords
76
 * This plugin uses the following keywords:
77
78
79
 * - COORDSYS: coordinate system (see Coords)
 * - COORDPARAM: extra coordinates parameters (see Coords)
 * - GRID2DPATH: path where to find the grids
80
 * - 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
81
 * - GRID2DEXT: grib file extension, or <i>none</i> for no file extension (default: ".grb")
82
83
 * - 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)
84
 * - METEOPATH: path where to find the grids for extracting time series at special points
Mathias Bavay's avatar
Mathias Bavay committed
85
 * - METEOEXT: file extension, or <i>none</i> for no file extension (default: ".grb")
86
 * - STATION#: coordinates for virtual stations (if using GRIB as METEO plugin). Each station is given by its coordinates and the closest
87
88
 * 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.
89
 *
90
91
92
 */

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

96
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),
          cellsize_x(IOUtils::nodata), cellsize_y(IOUtils::nodata),
          indexed(false), meteo_initialized(false), update_dem(false)
103
104
105
106
{
	setOptions();
}

107
108
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),
          cellsize_x(IOUtils::nodata), cellsize_y(IOUtils::nodata),
          indexed(false), meteo_initialized(false), update_dem(false)
114
115
{
	setOptions();
116
117
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
143
144
145
}

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;
		indexed = source.indexed;
		meteo_initialized = source.meteo_initialized;
		update_dem = source.update_dem;
	}
	return *this;
146
147
148
149
}

GRIBIO::~GRIBIO() throw()
{
150
	cleanup();
151
152
153
154
}

void GRIBIO::setOptions()
{
155
	std::string coordout, coordoutparam;
156
	IOUtils::getProjectionParameters(cfg, coordin, coordinparam, coordout, coordoutparam);
157

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

166
	cfg.getValue("METEOEXT", "Input", meteo_ext, IOUtils::nothrow);
167
	if(meteo_ext=="none") meteo_ext.clear();
168

169
	cfg.getValue("GRID2DEXT", "Input", grid2d_ext, IOUtils::nothrow);
170
	if(grid2d_ext=="none") grid2d_ext.clear();
171
}
172

173
void GRIBIO::readStations(std::vector<Coords> &vecPoints)
174
175
176
177
178
{
	cfg.getValue("METEOPATH", "Input", meteopath_in);
	size_t current_stationnr = 1;
	string current_station;
	do {
179
		current_station.clear();
180
		ostringstream ss;
181
		ss << "STATION" << current_stationnr;
182
		cfg.getValue(ss.str(), "Input", current_station, IOUtils::nothrow);
183
		IOUtils::stripComments(current_station);
184

185
		if (!current_station.empty()) {
186
187
			addStation(current_station, vecPoints);
			std::cerr <<  "\tRead virtual station " << vecPoints.back().printLatLon() << "\n";
188
189
		}
		current_stationnr++;
190
	} while (!current_station.empty());
191
192
}

193
void GRIBIO::addStation(const std::string& coord_spec, std::vector<Coords> &vecPoints)
194
195
196
197
198
199
200
201
202
203
204
205
206
{
	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);
207
		vecPoints.push_back(point);
208
209
210
211
212
		return;
	}
	if(coord1!=IOUtils::nodata && coord2!=IOUtils::nodata) {
		Coords point(coordin, coordinparam);
		point.setLatLon(coord1, coord2, IOUtils::nodata);
213
		vecPoints.push_back(point);
214
215
216
217
		return;
	}

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

220
void GRIBIO::listKeys(grib_handle** h, const std::string& filename)
221
{
222
223
224
225
	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);
226
227
228
229
230

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

232
233
234
	//Iterating over all keys
	while(grib_keys_iterator_next(kiter)) {
		const char* name = grib_keys_iterator_get_name(kiter);
235
236
		char value[MAX_VAL_LEN];
		size_t vlen = MAX_VAL_LEN;
237
238
		bzero(value,vlen);
		GRIB_CHECK(grib_get_string(*h,name,value,&vlen),name);
239
		std::cerr << name << " = " << value << "\n";
240
241
242
243
244
	}

	grib_keys_iterator_delete(kiter);
}

245
void  GRIBIO::listFields(const std::string& filename)
246
{
247
	grib_handle* h=NULL;
248
	int err=0;
249

250
	//loop over the messages
251
252
253
254
255
256
	while((h = grib_handle_new_from_file(0,fp,&err)) != NULL) {
		if(!h) {
			cleanup();
			throw IOException("Unable to create grib handle for \""+filename+"\"", AT);
		}

257
258
259
260
261
262
263
264
		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);
265
		long levelType;
266
267
268
		GRIB_CHECK(grib_get_long(h,"indicatorOfTypeOfLevel", &levelType),0);
		long level=0;
		if(levelType!=1) GRIB_CHECK(grib_get_long(h,"level", &level),0);
269
		std::cerr << marsParam << " " << shortname << " " << name << " type " << levelType << " level " << level << "\n";
270
		grib_handle_delete(h);
271
272
	}
}
273

274
void GRIBIO::getDate(grib_handle* h, Date &base, double &d1, double &d2) {
275
276
277
278
	long dataDate, dataTime;
	GRIB_CHECK(grib_get_long(h,"dataDate",&dataDate),0);
	GRIB_CHECK(grib_get_long(h,"dataTime",&dataTime),0);

279
280
	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);
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
307
308
309
310
311
312
313
314
315
316
317
318
	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:
319
			std::ostringstream ss;
320
321
322
323
			ss << "GRIB file using stepUnits=" << stepUnits << ", which is not supported";
			throw InvalidFormatException(ss.str(), AT);
	}

324
325
	d1 = static_cast<double>(startStep*step_units);
	d2 = static_cast<double>(endStep*step_units);
326
327
}

328
Coords GRIBIO::getGeolocalization(grib_handle* h, double &cell_x, double &cell_y)
329
{
330
	//getting transformation parameters
331
332
333
334
	double angleOfRotationInDegrees;
	GRIB_CHECK(grib_get_double(h,"angleOfRotationInDegrees",&angleOfRotationInDegrees),0);
	if(angleOfRotationInDegrees!=0.) {
		throw InvalidArgumentException("Rotated grids not supported!", AT);
335
	}
336
	double latitudeOfSouthernPole, longitudeOfSouthernPole;
337
338
	GRIB_CHECK(grib_get_double(h,"latitudeOfSouthernPoleInDegrees",&latitudeOfSouthernPole),0);
	GRIB_CHECK(grib_get_double(h,"longitudeOfSouthernPoleInDegrees",&longitudeOfSouthernPole),0);
339
340
	latitudeOfNorthernPole = -latitudeOfSouthernPole;
	longitudeOfNorthernPole = longitudeOfSouthernPole+180.;
341

342
343
	//determining llcorner, urcorner and center coordinates
	double ll_latitude, ll_longitude, ll_lat, ll_lon;
344
345
	GRIB_CHECK(grib_get_double(h,"latitudeOfFirstGridPointInDegrees",&ll_latitude),0);
	GRIB_CHECK(grib_get_double(h,"longitudeOfFirstGridPointInDegrees",&ll_longitude),0);
346
347
	Coords::rotatedToTrueLatLon(latitudeOfNorthernPole, longitudeOfNorthernPole, ll_latitude, ll_longitude, ll_lat, ll_lon);
	double ur_latitude, ur_longitude, ur_lat, ur_lon;
348
349
	GRIB_CHECK(grib_get_double(h,"latitudeOfLastGridPointInDegrees",&ur_latitude),0);
	GRIB_CHECK(grib_get_double(h,"longitudeOfLastGridPointInDegrees",&ur_longitude),0);
350
	Coords::rotatedToTrueLatLon(latitudeOfNorthernPole, longitudeOfNorthernPole, ur_latitude, ur_longitude, ur_lat, ur_lon);
351
	double cntr_lat, cntr_lon; //geographic coordinates
352
	Coords::rotatedToTrueLatLon(latitudeOfNorthernPole, longitudeOfNorthernPole, .5*(ll_latitude+ur_latitude), .5*(ll_longitude+ur_longitude), cntr_lat, cntr_lon);
353
354

	//determining cell size
355
356
357
	long Ni, Nj;
	GRIB_CHECK(grib_get_long(h,"Ni",&Ni),0);
	GRIB_CHECK(grib_get_long(h,"Nj",&Nj),0);
358
	double bearing;
359
360
361
362
363
364
365
366
	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)
367

368
369
370
	//computing lower left corner by using the center point as reference
	Coords cntr(coordin, coordinparam);
	cntr.setLatLon(cntr_lat, cntr_lon, IOUtils::nodata);
371
372
	const double cellsize=.5*(cell_x+cell_y);
	cntr.moveByXY(-.5*(double)Ni*cellsize, -.5*(double)Nj*cellsize);
373

374
	//checking that cellsize does not vary too much across the grid
375
	/*const double cellsize_x_ll = Coords::lon_degree_lenght(ll_latitude)*d_j;
376
377
	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.) {
378
		ostringstream ss;
379
380
381
		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);
382
	}*/
383

384
	return cntr; //this is now the ll corner
385
386
}

387
void GRIBIO::read2Dlevel(grib_handle* h, Grid2DObject& grid_out, const bool& read_geolocalization)
388
{
389
390
391
392
393
394
395
	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)) {
396
		ostringstream ss;
397
398
399
400
		ss << "Declaring grid of size " << Ni << "x" << Nj << "=" << Ni*Nj << " ";
		ss << "but containing " << values_len << " values. This is inconsistent!";
		throw InvalidArgumentException(ss.str(), AT);
	}
401
	double *values = (double*)calloc(values_len, sizeof(double));
402
403

	GRIB_CHECK(grib_get_double_array(h,"values",values,&values_len),0);
404
405
406
407

	if(read_geolocalization) {
		llcorner = getGeolocalization(h, cellsize_x, cellsize_y);
		if( fabs(cellsize_x-cellsize_y)/cellsize_x > 1./100.) {
408
409
			free(values);
			throw InvalidArgumentException("Unsupported geometry: cells can not be represented by square cells!", AT);
410
		}
411
	}
412
413
414
415
	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++)
416
			grid_out(ii,jj) = values[i++];
417
	}
418
419
	free(values);
}
420

421
bool GRIBIO::read2DGrid_indexed(const double& in_marsParam, const long& i_levelType, const long& i_level, const Date i_date, Grid2DObject& grid_out)
422
{
423
424
425
	GRIB_CHECK(grib_index_select_double(idx,"marsParam",in_marsParam),0);
	GRIB_CHECK(grib_index_select_long(idx,"indicatorOfTypeOfLevel", i_levelType),0);

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

434
435
436
437
438
439
440
441
		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);

442
		long level=0;
443
		if(i_level!=0) GRIB_CHECK(grib_get_long(h,"level", &level),0);
444
445
446
447
448
449
450
451
452
453
		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;
			}
454
		}
455
		grib_handle_delete(h);
456
	}
457
	return false;
458
459
}

460
void GRIBIO::read2DGrid(Grid2DObject& grid_out, const std::string& i_name)
461
{
462
463
464
	const std::string filename = grid2dpath_in+"/"+i_name;
	fp = fopen(filename.c_str(),"r");
	if(fp==NULL) {
465
		ostringstream ss;
466
467
468
469
470
471
472
473
474
475
476
477
		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);
		}

478
		read2Dlevel(h, grid_out, true);
479
480
481
482
483
484
485
		grib_handle_delete(h);
	} else {
		cleanup();
		throw IOException("No grid found in file \""+filename+"\"", AT);
	}

	cleanup();
486
487
}

488
void GRIBIO::indexFile(const std::string& filename)
489
{
490
491
	fp = fopen(filename.c_str(),"r");
	if(fp==NULL) {
492
		ostringstream ss;
493
494
495
496
		ss << "Error openning file \"" << filename << "\", possible reason: " << strerror(errno);
		throw FileAccessException(ss.str(), AT);
	}

497
	int err=0;
498
	const std::string keys("marsParam:d,indicatorOfTypeOfLevel:l"); //indexing keys
499
500
501
502
	char *c_filename = (char *)filename.c_str();
	idx = grib_index_new_from_file(0, c_filename, keys.c_str(), &err);
	if(err!=0) {
		cleanup();
503
		throw IOException("Failed to index GRIB file \""+filename+"\". Is it a valid GRIB file?", AT);
504
505
506
	}
	indexed=true;
	idx_filename=filename;
507
508
509
510
511
512
513
514
515
516

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

517
		//llcorner = getGeolocalization(h, cellsize_x, cellsize_y); //this sets llcorner, cellsize and bearing_offset
518
519
520
521
522
523
524
525
526
		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);
	}

527
}
528

529
void GRIBIO::read2DGrid(Grid2DObject& grid_out, const MeteoGrids::Parameters& parameter, const Date& date)
530
{
531
532
	Date UTC_date = date;
	UTC_date.setTimeZone(tz_in);
533

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

	read2DGrid(filename, grid_out, parameter, UTC_date);
537
538
}

539
540
541
542
543
544
545
546
547
548
549
550
551
552
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);
553
554
		for(size_t jj=0; jj<VW.nrows; jj++) {
			for(size_t ii=0; ii<VW.ncols; ii++) {
555
				VW(ii,jj) = sqrt( Optim::pow2(U(ii,jj)) + Optim::pow2(V(ii,jj)) );
556
				DW(ii,jj) = fmod( atan2( U(ii,jj), V(ii,jj) ) * Cst::to_deg + 360. + bearing_offset, 360.); // turn into degrees [0;360)
557
558
559
560
561
562
563
			}
		}
	}

	wind_date = date;
}

564
void GRIBIO::read2DGrid(const std::string& filename, Grid2DObject& grid_out, const MeteoGrids::Parameters& parameter, const Date& date)
565
{ //Parameters should be read in table 2 if available since this table is the one standardized by WMO
566
	if(!indexed || idx_filename!=filename) {
567
568
569
570
		cleanup();
		indexFile(filename);
	}

571
	//Basic meteo parameters
572
	if(parameter==MeteoGrids::P) read2DGrid_indexed(1.2, 1, 0, date, grid_out); //PS
573
	if(parameter==MeteoGrids::TA) read2DGrid_indexed(11.2, 105, 2, date, grid_out); //T_2M
574
	if(parameter==MeteoGrids::RH) {
575
		if(!read2DGrid_indexed(52.2, 105, 2, date, grid_out)) { //RELHUM_2M
576
			Grid2DObject ta;
577
578
			read2DGrid_indexed(11.2, 105, 2, date, ta); //T_2M
			read2DGrid_indexed(17.2, 105, 2, date, grid_out); //TD_2M
579
580
			for(size_t jj=0; jj<grid_out.nrows; jj++) {
				for(size_t ii=0; ii<grid_out.ncols; ii++) {
581
582
583
584
585
					grid_out(ii,jj) = Atmosphere::DewPointtoRh(grid_out(ii,jj), ta(ii,jj), true);
				}
			}
		}
	}
586
587
	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
588
589

	//hydrological parameters
590
591
592
593
	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) {
594
595
596
597
598
599
		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;
		}
600
	}
601
602

	//radiation parameters
603
604
605
606
	if(parameter==MeteoGrids::ALB) {
		read2DGrid_indexed(84.2, 1, 0, date, grid_out); //ALB_RAD
		grid_out.grid2D /= 100.;
	}
607
608
609
610
611
	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
	}
612
	/*if(parameter==MeteoGrids::CLD) { //cloudiness
613
614
615
616
		if(read2DGrid_indexed(74.2, 1, 0, date, grid_out)) //CLCM
		grid_out.grid2D /= 100.;
	}*/

617
	if(parameter==MeteoGrids::ISWR) {
618
619
		if(read2DGrid_indexed(116.2, 1, 0, date, grid_out)) { //short wave
			grid_out.grid2D *= -1.;
620
		} else if(!read2DGrid_indexed(111.250, 1, 0, date, grid_out)) { //GLOB
621
			Grid2DObject diff;
622
623
			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
624
625
626
627
628
			grid_out.grid2D += diff.grid2D;
		}
	}

	//DEM parameters
629
	if(parameter==MeteoGrids::DEM) read2DGrid_indexed(8.2, 1, 0, date, grid_out); //HSURF
630
631
	if(parameter==MeteoGrids::SLOPE) {
		read2DGrid_indexed(98.202, 1, 0, date, grid_out); //SLO_ANG
632
		grid_out.grid2D *= Cst::to_deg;
633
634
635
	}
	if(parameter==MeteoGrids::AZI) {
		read2DGrid_indexed(99.202, 1, 0, date, grid_out); //SLO_ASP
636
637
		for(size_t jj=0; jj<grid_out.nrows; jj++) {
			for(size_t ii=0; ii<grid_out.ncols; ii++) {
638
				grid_out(ii,jj) = fmod( grid_out(ii,jj)*Cst::to_deg + 360. + bearing_offset, 360.); // turn into degrees [0;360)
639
640
			}
		}
641
	}
642

643
	//Wind parameters
644
	if(parameter==MeteoGrids::VW_MAX) read2DGrid_indexed(187.201, 105, 10, date, grid_out); //VMAX_10M 10m
645
646
647
648
	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);
649
650
		for(size_t jj=0; jj<grid_out.nrows; jj++) {
			for(size_t ii=0; ii<grid_out.ncols; ii++) {
651
				grid_out(ii,jj) = VW(ii,jj)*sin(DW(ii,jj)*Cst::to_rad);
652
653
654
			}
		}
	}
655
656
	if(parameter==MeteoGrids::V) {
		readWind(filename, date);
657
658
		for(size_t jj=0; jj<grid_out.nrows; jj++) {
			for(size_t ii=0; ii<grid_out.ncols; ii++) {
659
				grid_out(ii,jj) = VW(ii,jj)*cos(DW(ii,jj)*Cst::to_rad);
660
661
662
			}
		}
	}
663
664
665
666
667
668
669
670
	if(parameter==MeteoGrids::DW) {
		readWind(filename, date);
		grid_out = DW;
	}
	if(parameter==MeteoGrids::VW) {
		readWind(filename, date);
		grid_out = VW;
	}
671
672

	if(grid_out.isEmpty()) {
673
		ostringstream ss;
674
675
676
677
678
679
		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
680
681
	/*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
682
		Grid2DObject Z0;
683
		if(read2DGrid_indexed(83.2, 1, 0, date, Z0)) { //Z0
684
685
			for(size_t jj=0; jj<grid_out.nrows; jj++) {
				for(size_t ii=0; ii<grid_out.ncols; ii++) {
686
					grid_out(ii,jj) = Atmosphere::windLogProfile(grid_out(ii,jj), 10., 7.5, Z0(ii,jj));
687
688
				}
			}
689
		} else {
690
			const double wind_factor = Atmosphere::windLogProfile(1., 10., 7.5, 0.03);
691
			grid_out.grid2D *= wind_factor;
692
		}
693
	}*/
694
695
}

696
void GRIBIO::readDEM(DEMObject& dem_out)
697
{
698
	const Date d(2012,9,19,0,0,0); //ie: undef. This will be caught when reading the GRIB file
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
	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();
	}
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
}

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

747
748
749
void GRIBIO::scanMeteoPath()
{
	std::list<std::string> dirlist;
750
	IOUtils::readDirectory(meteopath_in, dirlist, meteo_ext);
751
752
753
754
755
756
757
758
	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;
759
		IOUtils::convertString(date, filename.substr(spos,10), tz_in);
760
761
762
763
764
765
766
		std::pair<Date,std::string> tmp(date, filename);

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

767
768
void GRIBIO::readMeteoData(const Date& dateStart, const Date& dateEnd,
                             std::vector< std::vector<MeteoData> >& vecMeteo,
769
770
                             const size_t&)
{
771
	if(!meteo_initialized) {
772
		readStations(vecPts);
773
774
775
		scanMeteoPath();
		meteo_initialized=true;
	}
776

777
	vecMeteo.clear();
778

779
780
	double *lats = (double*)malloc(vecPts.size()*sizeof(double));
	double *lons = (double*)malloc(vecPts.size()*sizeof(double));
781
782
783
784
	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
785
	size_t idx_start;
786
787
788
789
790
791
792
793
794
	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)
795
796

	try {
797
		for(size_t ii=idx_start; ii<cache_meteo_files.size(); ii++) {
798
799
800
			const Date& date = cache_meteo_files[ii].first;
			if(date>dateEnd) break;
			const std::string filename = meteopath_in+"/"+cache_meteo_files[ii].second;
801
802
803

			if(!indexed || idx_filename!=filename) {
				cleanup();
804
				indexFile(filename); //this will also read geolocalization
805
806
			}
			if(!meta_ok) {
807
				if(readMeteoMeta(vecPts, meta, lats, lons)==false) {
808
809
810
811
					//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));
812
					readMeteoMeta(vecPts, meta, lats, lons);
813
				}
814
				vecMeteo.insert(vecMeteo.begin(), vecPts.size(), std::vector<MeteoData>()); //allocation for the vectors now that we know how many true stations we have
815
				meta_ok=true;
816
817
818
819
			}

			std::vector<MeteoData> Meteo;
			readMeteoStep(meta, lats, lons, date, Meteo);
820
			for(size_t jj=0; jj<vecPts.size(); jj++)
821
822
823
824
				vecMeteo[jj].push_back(Meteo[jj]);
		}
	} catch(...) {
		free(lats); free(lons);
825
		cleanup();
826
		throw;
827
828
	}

829
	free(lats); free(lons);
830
831
}

832
bool GRIBIO::removeDuplicatePoints(std::vector<Coords>& vecPoints, double *lats, double *lons)
833
{ //remove potential duplicates. Returns true if some have been removed
834
	const size_t npoints = vecPoints.size();
835
836
	std::vector<size_t> deletions;
	deletions.reserve(npoints);
837
	for(size_t ii=0; ii<npoints; ii++) {
838
839
		const double lat = lats[ii];
		const double lon = lons[ii];
840
		for(size_t jj=ii+1; jj<npoints; jj++) {
841
842
843
			if(lat==lats[jj] && lon==lons[jj]) {
				deletions.push_back(jj);
			}
844
845
846
847
		}
	}

	//we need to erase from the end in order to keep the index unchanged...
848
849
	for(size_t ii=deletions.size(); ii>0; ii--) {
		const size_t index=deletions[ii-1];
850
		vecPoints.erase(vecPoints.begin()+index);
851
852
	}

853
	if(!deletions.empty()) return true;
854
855
856
	return false;
}

857
858
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)
859
860
	stations.clear();

861
862
863
864
	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;
865
866
867
868
869
	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);
	}
870

871
	const size_t npoints = vecPoints.size();
872
	double latitudeOfSouthernPole, longitudeOfSouthernPole;
873
874
	GRIB_CHECK(grib_get_double(h,"latitudeOfSouthernPoleInDegrees",&latitudeOfSouthernPole),0);
	GRIB_CHECK(grib_get_double(h,"longitudeOfSouthernPoleInDegrees",&longitudeOfSouthernPole),0);
875
876
	latitudeOfNorthernPole = -latitudeOfSouthernPole;
	longitudeOfNorthernPole = longitudeOfSouthernPole+180.;
877

878
879
	long Ni;
	GRIB_CHECK(grib_get_long(h,"Ni",&Ni),0);
880

881
	//build GRIB local coordinates for the points
882
	for(size_t ii=0; ii<npoints; ii++) {
883
		Coords::trueLatLonToRotated(latitudeOfNorthernPole, longitudeOfNorthernPole, vecPoints[ii].getLat(), vecPoints[ii].getLon(), lats[ii], lons[ii]);
884
	}
885

886
887
888
889
890
891
892
893
894
895
896
	//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);
	}
897

898
	//remove potential duplicates
899
	if(removeDuplicatePoints(vecPoints, outlats, outlons)==true) {
900
901
		free(outlats); free(outlons); free(values); free(distances); free(indexes);
		grib_handle_delete(h);
902
		return false;
903
	}
904
905

	//fill metadata
906
	for(size_t ii=0; ii<npoints; ii++) {
907
908
		StationData sd;
		sd.position.setProj(coordin, coordinparam);
909
		double true_lat, true_lon;
910
		Coords::rotatedToTrueLatLon(latitudeOfNorthernPole, longitudeOfNorthernPole, outlats[ii], outlons[ii], true_lat, true_lon);
911
		sd.position.setLatLon(true_lat, true_lon, values[ii]);
912
		ostringstream ss;
913
914
		ss << "Point_" << indexes[ii];
		sd.stationID=ss.str();
915
		ostringstream ss2;
916
917
918
919
920
921
922
923
		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;
924
925
}

926
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)
927
928
929
{
	GRIB_CHECK(grib_index_select_double(idx,"marsParam",marsParam),0);
	GRIB_CHECK(grib_index_select_long(idx,"indicatorOfTypeOfLevel", levelType),0);
930
931
932
933
934
935
936
937
938

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

939
940
941
942
943
944
945
946
		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);

947
948
		long level=0;
		if(i_level!=0) GRIB_CHECK(grib_get_long(h,"level", &level),0);
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
		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);
966
				grib_handle_delete(h);
967
				return true;
968
969
970
971
972
			}
		}
		grib_handle_delete(h);
	}
	return false;