WSL/SLF GitLab Repository

NetCDFIO.cc 36.4 KB
Newer Older
1
/***********************************************************************************/
2
/*  Copyright 2014 WSL Institute for Snow and Avalanche Research    SLF-DAVOS      */
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
/***********************************************************************************/
/* This file is part of MeteoIO.
    MeteoIO is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    MeteoIO is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public License
    along with MeteoIO.  If not, see <http://www.gnu.org/licenses/>.
*/
#include "NetCDFIO.h"
19
#include <meteoio/ResamplingAlgorithms2D.h>
20
#include <meteoio/meteoStats/libinterpol1D.h>
21
22
#include <meteoio/meteoLaws/Meteoconst.h>
#include <meteoio/meteoLaws/Atmosphere.h>
23
#include <meteoio/Timer.h>
24
#include <meteoio/MathOptim.h>
25
26
27
28
29
#include <meteoio/plugins/libncpp.h>

#include <cmath>
#include <cstdio>
#include <algorithm>
30
31
32
33
34
35
36

using namespace std;

namespace mio {
/**
 * @page netcdf NetCDF
 * @section netcdf_format Format
37
38
39
40
41
 * In order to promote creation, access and sharing of scientific data, the NetCDF format has been
 * created as a machine-independent format. NetCDF (network Common Data Form) is therefore an interface
 * for array-oriented data access and a library that provides an implementation of the interface. The
 * <A HREF="http://www.unidata.ucar.edu/downloads/netcdf/index.jsp">NetCDF software</A> was developed
 * at the <A HREF="http://www.unidata.ucar.edu/">Unidata Program Center</A> in Boulder, Colorado.
42
 * In order to graphicaly explore the content and structure of NetCDF files, you can use the
43
 * <A HREF="http://www.epic.noaa.gov/java/ncBrowse/">ncBrowse</A> java software.
44
 *
45
46
47
48
 * The NetCDF format does not impose a specific set of metadata and therefore in order to easily exchange data
 * within a given field, it is a good idea to standardize the metadata. Several such metadata schema can be used
 * by this plugin:
 * - CF1 - the <A HREF="http://cfconventions.org">conventions</A> for climate and forecast (CF) metadata;
Mathias Bavay's avatar
Mathias Bavay committed
49
 * - ECMWF - from the <A HREF="http://www.ecmwf.int/">European Centre for Medium-Range Weather Forecasts</A>, see the <A HREF="https://software.ecmwf.int/wiki/display/TIGGE/Soil+temperature">ECMWF Wiki</A> for a description of the available fields;
50
 * - CNRM - from the <A HREF="http://www.cnrm.meteo.fr/">National Centre for Meteorological Research</A>.
51
52
53
54
55
56
57
58
 *
 * @section netcdf_units Units
 *
 *
 * @section netcdf_keywords Keywords
 * This plugin uses the following keywords:
 * - COORDSYS: coordinate system (see Coords); [Input] and [Output] section
 * - COORDPARAM: extra coordinates parameters (see Coords); [Input] and [Output] section
59
60
 * - DEMFILE: The filename of the file containing the DEM; [Input] section
 * - DEMVAR: The variable name of the DEM within the DEMFILE; [Input] section
Mathias Bavay's avatar
Mathias Bavay committed
61
62
63
64
 * - GRID2DPATH: if this directory contains files, they will be used for reading the input from; [Input] section
 * - METEO_EXT: only the files containing this pattern in their filename will be used; [Input] section
 * - GRID2DFILE: if GRID2DPATH has not been defined or if it does not contain files matching the METEO_EXT extension, provides
 * the NetCDF file which shall be used for gridded input/output; [Input] and [Output] section
65
 * - NETCDF_SCHEMA: the schema to use (either CF1 or CNRM or ECMWF); [Input] and [Output] section
66
 * - NETCDF::{MeteoGrids::Parameters} = {netcdf_param_name} : this allows to remap the names as found in the NetCDF file to the MeteoIO grid parameters; [Input] section;
67
 * - DEM_FROM_PRESSURE: if no dem is found but local and sea level pressure grids are found, use them to rebuild a DEM; [Input] section
68
 *
Mathias Bavay's avatar
Mathias Bavay committed
69
70
71
 * When providing multiple files in one directory, in case of overlapping files, the file containing the newest data has priority. This is
 * convenient when using forecats data to automatically use the most short-term forecast.
 * 
72
 * @section netcdf_example Example use
73
74
75
76
 * @code
 * [Input]
 * DEM     = NETCDF
 * DEMFILE = ./input/Aster_tile.nc
77
 * ;DEMVAR = height        ;this is only required if the variable name is non-standard
Mathias Bavay's avatar
Mathias Bavay committed
78
79
80
81
82
 * 
 * GRID2D    = NETCDF
 * GRID2DPATH =  /data/meteo_reanalysis
 * METEO_EXT = .nc
 * NETCDF_SCHEMA = ECMWF
83
 * NETCDF::PSUM = RhiresD               ;overwrite the PSUM parameter with "RhiresD", for example for MeteoCH reanalysis
84
85
 * @endcode
 *
86
 * @section netcdf_compilation Compilation
87
88
 * In order to compile this plugin, you need libnetcdf (for C). For Linux, please select both the libraries and
 * their development files in your package manager.
89
90
 */

91
const double NetCDFIO::plugin_nodata = -9999999.; //CNRM-GAME nodata value
92
const std::string NetCDFIO::cf_time = "time";
93
94
const std::string NetCDFIO::cf_latitude = "latitude";
const std::string NetCDFIO::cf_longitude = "longitude";
95

96
97
NetCDFIO::NetCDFIO(const std::string& configfile) : cfg(configfile), cache_meteo_files(), in_attributes(), out_attributes(), 
                                                    coordin(), coordinparam(), coordout(), coordoutparam(), 
98
                                                    in_dflt_TZ(0.), out_dflt_TZ(0.), in_time_offset(IOUtils::nodata), in_time_multiplier(IOUtils::nodata), 
99
                                                    dem_altimeter(false), in_strict(false), out_strict(false), vecMetaData()
100
101
{
	IOUtils::getProjectionParameters(cfg, coordin, coordinparam, coordout, coordoutparam);
102
	parseInputOutputSection();
103
104
}

105
106
NetCDFIO::NetCDFIO(const Config& cfgreader) : cfg(cfgreader), cache_meteo_files(), in_attributes(), out_attributes(), 
                                              coordin(), coordinparam(), coordout(), coordoutparam(), 
107
                                              in_dflt_TZ(0.), out_dflt_TZ(0.), in_time_offset(IOUtils::nodata), in_time_multiplier(IOUtils::nodata), 
108
                                              dem_altimeter(false), in_strict(false), out_strict(false), vecMetaData()
109
110
{
	IOUtils::getProjectionParameters(cfg, coordin, coordinparam, coordout, coordoutparam);
111
	parseInputOutputSection();
112
113
}

114
NetCDFIO::~NetCDFIO() throw() {}
115

116
117
118
119
120
121
void NetCDFIO::parseInputOutputSection()
{
	//default timezones
	in_dflt_TZ = out_dflt_TZ = IOUtils::nodata;
	cfg.getValue("TIME_ZONE", "Input", in_dflt_TZ, IOUtils::nothrow);
	cfg.getValue("TIME_ZONE", "Output", out_dflt_TZ, IOUtils::nothrow);
122
	cfg.getValue("DEM_FROM_PRESSURE", "Input", dem_altimeter, IOUtils::nothrow);
123
	
124
125
	string in_schema = "ECMWF";
	in_schema = IOUtils::strToUpper( cfg.get("NETCDF_SCHEMA", "Input", IOUtils::nothrow) );
126
	initAttributesMap(in_schema, in_attributes);
127
128
	string out_schema = "ECMWF";
	out_schema = IOUtils::strToUpper( cfg.get("NETCDF_SCHEMA", "Output", IOUtils::nothrow) );
129
	initAttributesMap(out_schema, out_attributes);
130
131
132
133
	
	string in_grid2d_path;
	cfg.getValue("GRID2DPATH", "Input", in_grid2d_path, IOUtils::nothrow);
	if (!in_grid2d_path.empty()) scanMeteoPath(in_grid2d_path,  cache_meteo_files);
134
135
}

136
void NetCDFIO::initAttributesMap(const std::string& schema, std::map<MeteoGrids::Parameters, attributes> &attr)
137
138
139
140
141
142
143
144
145
146
147
148
149
150
{
	if (schema.empty()) return;
	
	if (schema=="CF1") {
		attr[MeteoGrids::DEM] = attributes("z", "altitude", "height above mean sea level", "m", IOUtils::nodata);
		attr[MeteoGrids::TA] = attributes("temperature", "air_temperature", "near surface air temperature", "K", IOUtils::nodata);
		attr[MeteoGrids::RH] = attributes("humidity", "relative humidity", "relative humidity", "fraction", IOUtils::nodata);
		attr[MeteoGrids::P] = attributes("pressure", "air_pressure", "near surface air pressure", "Pa", IOUtils::nodata);
	} else if (schema=="CNRM") {
		attr[MeteoGrids::DEM] = attributes("ZS", "", "altitude", "m", IOUtils::nodata);
		attr[MeteoGrids::SLOPE] = attributes("slope", "", "slope angle", "degrees from horizontal", IOUtils::nodata);
		attr[MeteoGrids::AZI] = attributes("aspect", "", "slope aspect", "degrees from north", IOUtils::nodata);
		attr[MeteoGrids::TA] = attributes("Tair", "", "Near Surface Air Temperature", "K", IOUtils::nodata);
		attr[MeteoGrids::RH] = attributes("HUMREL", "", "Relative Humidity", "%", IOUtils::nodata);
151
		attr[MeteoGrids::QI] = attributes("Qair", "", "", "", IOUtils::nodata);
152
153
		attr[MeteoGrids::VW] = attributes("Wind", "", "Wind Speed", "m/s", IOUtils::nodata);
		attr[MeteoGrids::DW] = attributes("Wind_DIR", "", "Wind Direction", "deg", IOUtils::nodata);
154
155
		attr[MeteoGrids::PSUM_L] = attributes("Rainf", "", "Rainfall Rate", "kg/m2/s", IOUtils::nodata);
		attr[MeteoGrids::PSUM_S] = attributes("Snowf", "", "", "", IOUtils::nodata);
156
157
		attr[MeteoGrids::ISWR_DIR] = attributes("DIR_SWdown", "", "Surface Incident Direct Shortwave Radiation", "W/m2", IOUtils::nodata);
		attr[MeteoGrids::ISWR_DIFF] = attributes("SCA_SWdown", "", "", "", IOUtils::nodata);
158
159
160
		attr[MeteoGrids::P] = attributes("PSurf", "", "Surface Pressure", "Pa", IOUtils::nodata);
		attr[MeteoGrids::ILWR] = attributes("LWdown", "", "Surface Incident Longwave Radiation", "W/m2", IOUtils::nodata);
	} else if (schema=="ECMWF") {
161
		attr[MeteoGrids::DEM] = attributes("z", "geopotential_height", "geopotential_height", "m", IOUtils::nodata);
162
163
164
165
166
167
		attr[MeteoGrids::TA] = attributes("t2m", "", "2 metre temperature", "K", 2.);
		attr[MeteoGrids::TD] = attributes("d2m", "", "2 metre dewpoint temperature", "K", 2.);
		attr[MeteoGrids::P] = attributes("sp", "surface_air_pressure", "Surface pressure", "Pa", IOUtils::nodata);
		attr[MeteoGrids::P_SEA] = attributes("msl", "air_pressure_at_sea_level", "Mean sea level pressure", "Pa", IOUtils::nodata);
		attr[MeteoGrids::ISWR] = attributes("ssrd", "surface_downwelling_shortwave_flux_in_air", "Surface solar radiation downwards", "J m**-2", IOUtils::nodata);
		attr[MeteoGrids::ILWR] = attributes("strd", "", "Surface thermal radiation downwards", "J m**-2", IOUtils::nodata);
168
		attr[MeteoGrids::PSUM] = attributes("tp", "", "Total precipitation", "m", IOUtils::nodata);
169
170
171
172
		attr[MeteoGrids::U] = attributes("u10", "", "10 metre U wind component", "m s**-1", 10.);
		attr[MeteoGrids::V] = attributes("v10", "", "10 metre V wind component", "m s**-1", 10.);
		attr[MeteoGrids::SWE] = attributes("sd", "lwe_thickness_of_surface_snow_amount", "Snow depth", "m of water equivalent", IOUtils::nodata);
		attr[MeteoGrids::TSS] = attributes("skt", "", "Skin temperature", "K", IOUtils::nodata);
173
174
175
176
		attr[MeteoGrids::TSG] = attributes("stl1", "surface_temperature", "Soil temperature level 1", "K", IOUtils::nodata); //this is from 0 to -7cm
		attr[MeteoGrids::ALB] = attributes("al", "surface_albedo", "Albedo", "(0 - 1)", IOUtils::nodata);
		attr[MeteoGrids::RSNO] = attributes("rsn", "", "Snow density", "kg m**-3", IOUtils::nodata);
		attr[MeteoGrids::ROT] = attributes("ro", "", "Runoff", "m", IOUtils::nodata);
177
178
	} else
		throw InvalidArgumentException("Invalid schema selected for NetCDF: \""+schema+"\"", AT);
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
	
	vector<string> custom_attr;
	const size_t nrOfCustoms = cfg.findKeys(custom_attr, "NETCDF::", "Input");
	for (size_t ii=0; ii<nrOfCustoms; ++ii) {
		const size_t found = custom_attr[ii].find_last_of(":");
		if (found==std::string::npos || found==custom_attr[ii].length()) continue;

		const string meteo_grid = custom_attr[ii].substr(found+1);
		const string netcdf_param = cfg.get(custom_attr[ii], "Input");
		const size_t param_index = MeteoGrids::getParameterIndex(meteo_grid);
		if (param_index==IOUtils::npos)
			throw InvalidArgumentException("Parameter '"+meteo_grid+"' is not a valid MeteoGrid! Please correct key '"+custom_attr[ii]+"'", AT);
		
		attr[ static_cast<MeteoGrids::Parameters>(param_index) ] = attributes(netcdf_param, "", "", "", IOUtils::nodata);
	}
194
195
}

196
void NetCDFIO::read2DGrid(Grid2DObject& grid_out, const std::string& arguments)
197
{
198
199
200
201
202
203
204
205
	vector<string> vec_argument;
	IOUtils::readLineToVec(arguments, vec_argument, ':');

	if (vec_argument.size() == 2) {
		read2DGrid_internal(grid_out, vec_argument[0], vec_argument[1]);
	} else {
		throw InvalidArgumentException("The format for the arguments to NetCDFIO::read2DGrid is filename:varname", AT);
	}
206
207
}

208
void NetCDFIO::read2DGrid(Grid2DObject& grid_out, const MeteoGrids::Parameters& parameter, const Date& date)
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
{
	if (!cache_meteo_files.empty()) {
		for (size_t ii=0; ii<cache_meteo_files.size(); ii++) {
			const Date date_start = cache_meteo_files[ii].first.first;
			const Date date_end = cache_meteo_files[ii].first.second;
			if (date>=date_start && date<=date_end) {
				const string filename = cache_meteo_files[ii].second;
				read2DGrid(grid_out, parameter, date, filename);
				return;
			}
		}
		//the date was not found
		string in_grid2d_path;
		cfg.getValue("GRID2DPATH", "Input", in_grid2d_path);
		throw InvalidArgumentException("No Gridded data found for "+date.toString(Date::ISO)+"in '"+in_grid2d_path+"'", AT);
	} else {
		const string filename = cfg.get("GRID2DFILE", "Input");
		read2DGrid(grid_out, parameter, date, filename);
	}
}

void NetCDFIO::read2DGrid(Grid2DObject& grid_out, const MeteoGrids::Parameters& parameter, const Date& date, const std::string& filename)
231
{
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
	grid_out.clear();
	if (read2DGrid_internal(grid_out, filename, parameter, date)) return; //schema naming
	
	if (parameter==MeteoGrids::VW || parameter==MeteoGrids::DW) {	//VW, DW
		Grid2DObject U,V;
		const bool hasU = read2DGrid_internal(U, filename, MeteoGrids::U, date);
		const bool hasV = read2DGrid_internal(V, filename, MeteoGrids::V, date);
		if (hasU==true && hasV==true) {
			grid_out.set(U, IOUtils::nodata);
			if (parameter==MeteoGrids::VW) {
				for (size_t ii=0; ii<(grid_out.getNx()*grid_out.getNy()); ii++)
					grid_out(ii) = sqrt( Optim::pow2(U(ii)) + Optim::pow2(V(ii)) );
			} else {
				for (size_t ii=0; ii<(grid_out.getNx()*grid_out.getNy()); ii++)
					grid_out(ii) =  fmod( atan2( U(ii), V(ii) ) * Cst::to_deg + 360., 360.); // turn into degrees [0;360)
			}
			return;
		}
	}
	
	if (parameter==MeteoGrids::RH) {								//RH
		Grid2DObject ta;
		DEMObject dem;
		bool hasDEM = false;
		try {
			readDEM(dem);
			hasDEM = true;
		} catch(...){}
		const bool hasQI = read2DGrid_internal(grid_out, filename, MeteoGrids::QI, date);
		const bool hasTA = read2DGrid_internal(ta, filename, MeteoGrids::TA, date);
		if (hasQI && hasDEM && hasTA) {
263
			for (size_t ii=0; ii<(grid_out.getNx()*grid_out.getNy()); ii++)
264
265
266
267
268
269
				grid_out(ii) = Atmosphere::specToRelHumidity(dem(ii), ta(ii), grid_out(ii));
			return;
		}
		
		const bool hasTD = read2DGrid_internal(grid_out, filename, MeteoGrids::TD, date);
		if (hasTA && hasTD) {
270
			for (size_t ii=0; ii<(grid_out.getNx()*grid_out.getNy()); ii++)
271
272
273
274
275
				grid_out(ii) = Atmosphere::DewPointtoRh(grid_out(ii), ta(ii), false);
			return;
		}
	}
	
276
277
278
279
280
281
282
283
284
285
286
287
288
289
	if (parameter==MeteoGrids::RSWR) {								//RSWR
		bool hasISWR = false;
		try {
			read2DGrid(grid_out, MeteoGrids::ISWR, date);
			hasISWR = true;
		} catch(...){}
		Grid2DObject alb;
		const bool hasALB = read2DGrid_internal(alb, filename, MeteoGrids::ALB);
		if (hasALB && hasISWR) {
			grid_out *= alb;
			return;
		}
	}
	
290
291
292
293
294
295
296
297
298
299
	if (parameter==MeteoGrids::ISWR) {								//ISWR
		Grid2DObject iswr_diff;
		const bool hasISWR_DIFF = read2DGrid_internal(iswr_diff, filename, MeteoGrids::ISWR_DIFF);
		const bool hasISWR_DIR = read2DGrid_internal(grid_out, filename, MeteoGrids::ISWR_DIR);
		if (hasISWR_DIFF && hasISWR_DIR) {
			grid_out += iswr_diff;
			return;
		}
	}
	
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
	if (parameter==MeteoGrids::PSUM) {								//PSUM
		Grid2DObject psum_s;
		const bool hasPSUM_S = read2DGrid_internal(psum_s, filename, MeteoGrids::PSUM_S);
		const bool hasPSUM_L = read2DGrid_internal(grid_out, filename, MeteoGrids::PSUM_L);
		if (hasPSUM_S && hasPSUM_L) {
			grid_out += psum_s;
			return;
		}
	}
	
	if (parameter==MeteoGrids::PSUM_PH) {								//PSUM_PH
		Grid2DObject psum_l;
		const bool hasPSUM_S = read2DGrid_internal(grid_out, filename, MeteoGrids::PSUM_S);
		const bool hasPSUM_L = read2DGrid_internal(psum_l, filename, MeteoGrids::PSUM_L);
		if (hasPSUM_S && hasPSUM_L) {
			grid_out += psum_l;
			const size_t nrCells = grid_out.getNx()*grid_out.getNy();
			for (size_t ii=0; ii<nrCells; ii++) {
				const double psum = grid_out(ii);
				if (psum!=IOUtils::nodata && psum>0)
					grid_out(ii) = psum_l(ii) / psum;
			}
322
323
324
325
			return;
		}
	}
	
326
327
328
329
330
331
332
333
334
335
336
	if (parameter==MeteoGrids::HS) {								//HS
		Grid2DObject rsno;
		const bool hasRSNO = read2DGrid_internal(rsno, filename, MeteoGrids::RSNO);
		const bool hasSWE = read2DGrid_internal(grid_out, filename, MeteoGrids::SWE);
		if (hasRSNO && hasSWE) {
			grid_out *= 1000.0; //convert mm=kg/m^3 into kg
			grid_out /= rsno;
			return;
		}
	}
	
337
	throw InvalidArgumentException("Parameter \'"+MeteoGrids::getParameterName(parameter)+"\' either not found in file \'"+filename+"\' or not found in current NetCDF schema", AT);
338
339
}

340
void NetCDFIO::readDEM(DEMObject& dem_out)
341
{
342
343
344
345
346
347
	const string filename = cfg.get("DEMFILE", "Input");
	const string varname = cfg.get("DEMVAR", "Input", IOUtils::nothrow);
	if (!varname.empty()) {
		if (!read2DGrid_internal(dem_out, filename, varname))
			throw InvalidArgumentException("Variable \'"+varname+"\' not found in file \'"+filename+"\'", AT);
	} else {
348
		if (read2DGrid_internal(dem_out, filename, MeteoGrids::DEM)) return; //schema naming
349
350
351
		if (read2DGrid_internal(dem_out, filename, "Band1")) return; //ASTER naming
		if (read2DGrid_internal(dem_out, filename, "z")) return; //GDAL naming
		
352
		//last chance: read from pressure grids
353
354
		if (dem_altimeter) {
			Grid2DObject p, ta, p_sea;
355
356
357
			if (read2DGrid_internal(p_sea, filename, MeteoGrids::P_SEA) &&
			read2DGrid_internal(p, filename, MeteoGrids::P) && 
			read2DGrid_internal(ta, filename, MeteoGrids::TA)) {
358
				dem_out.set(p, IOUtils::nodata);
359
360
				const double k = Cst::gravity / (Cst::dry_adiabatique_lapse_rate * Cst::gaz_constant_dry_air);
				const double k_inv = 1./k;
361
				for (size_t ii=0; ii<(dem_out.getNx()*dem_out.getNy()); ii++) {
362
363
364
					const double K = pow(p(ii)/p_sea(ii), k_inv);
					dem_out(ii) = ta(ii)*Cst::earth_R0*(1.-K) / (Cst::dry_adiabatique_lapse_rate * Cst::earth_R0 - ta(ii)*(1.-K));
				}
365
366
				return;
			}
367
368
		}
		
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
		throw InvalidArgumentException("The variable containing the DEM could not be found. Please specify it using the DEMVAR key.", AT);
	}
}

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

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

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

void NetCDFIO::readMeteoData(const Date&, const Date&, std::vector< std::vector<MeteoData> >&, const size_t&)
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

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

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

void NetCDFIO::write2DGrid(const Grid2DObject& grid_in, const std::string& arguments)
{
	// arguments is a string of the format filname:varname
	vector<string> vec_argument;
413
	if (IOUtils::readLineToVec(arguments, vec_argument, ':')  != 2)
414
		throw InvalidArgumentException("The format for the arguments to NetCDFIO::write2DGrid is filename:varname", AT);
415

416
417
418
	const string name = vec_argument[1];
	const attributes attr(name, name, name, "", IOUtils::nodata);
	write2DGrid_internal(grid_in, vec_argument[0], attr);
419
420
421
422
423
}

void NetCDFIO::write2DGrid(const Grid2DObject& grid_in, const MeteoGrids::Parameters& parameter, const Date& date)
{
	const string filename = cfg.get("GRID2DFILE", "Output");
424
425
426
	
	const std::map<MeteoGrids::Parameters, attributes>::const_iterator it = in_attributes.find(parameter);
	if (it!=in_attributes.end()) {
427
		const bool isPrecip = (parameter==MeteoGrids::PSUM || parameter==MeteoGrids::SWE);
428
429
430
431
432
433
		write2DGrid_internal(grid_in, filename, it->second, date, isPrecip);
	} else {
		const string name = MeteoGrids::getParameterName(parameter);
		const attributes attr(name, name, name, "", IOUtils::nodata);
		write2DGrid_internal(grid_in, filename, attr, date);
	}
434
435
}

436
437
438
439
440
441
442
443
444
445
446
447
448
//custom function for sorting cache_meteo_files
struct sort_pred {
	bool operator()(const std::pair<std::pair<Date,Date>,string> &left, const std::pair<std::pair<Date,Date>,string> &right) {
		if (left.first.first < right.first.first) return true;
		if (left.first.first > right.first.first) return false;
		return left.first.second < right.first.second; //date_start equallity case
	}
};

void NetCDFIO::scanMeteoPath(const std::string& meteopath_in,  std::vector< std::pair<std::pair<mio::Date, mio::Date>, std::string> > &meteo_files)
{
	meteo_files.clear();

Mathias Bavay's avatar
Mathias Bavay committed
449
	const string meteo_ext = 	cfg.get("METEO_EXT", "INPUT", IOUtils::nothrow);
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
	std::list<std::string> dirlist = IOUtils::readDirectory(meteopath_in, meteo_ext);
	if (dirlist.empty()) return; //nothing to do if the directory is empty, we will transparently swap to using GRID2DFILE
	dirlist.sort();

	//Check date range in every filename and cache it
	std::list<std::string>::const_iterator it = dirlist.begin();
	while ((it != dirlist.end())) {
		const std::string& filename = meteopath_in + "/" + *it;
		
		if (!IOUtils::fileExists(filename)) throw AccessException(filename, AT); //prevent invalid filenames
		int ncid;
		ncpp::open_file(filename, NC_NOWRITE, ncid);
		getTimeTransform(ncid, in_time_offset, in_time_multiplier); //always re-read offset and multiplier, it might be different for each file
		double min, max;
		const bool status = ncpp::get_dimensionMinMax(ncid, NetCDFIO::cf_time, min, max);
		ncpp::close_file(filename, ncid); //no need to keep file open anymore
		
		if (!status) throw IOException("Could not get min/max time for file '"+filename+"'", AT);
		
		const Date d_min(min*in_time_multiplier + in_time_offset, in_dflt_TZ);
		const Date d_max(max*in_time_multiplier + in_time_offset, in_dflt_TZ);
		const std::pair<Date, Date> dateRange(d_min, d_max);
		const std::pair<std::pair<Date, Date>,std::string> tmp(dateRange, filename);
		meteo_files.push_back(tmp);
		it++;
	}
	std::sort(meteo_files.begin(), meteo_files.end(), sort_pred());
	
	//now handle overlaping files: truncate the end date of the file starting earlier
	for (size_t ii=0; ii<(meteo_files.size()-1); ii++) {
		if (meteo_files[ii].first.second > meteo_files[ii+1].first.first)
			meteo_files[ii].first.second = meteo_files[ii+1].first.first;
	}
}

485
486
487
488
489
bool NetCDFIO::read2DGrid_internal(Grid2DObject& grid_out, const std::string& filename, const MeteoGrids::Parameters& parameter, const Date& date)
{
	const std::map<MeteoGrids::Parameters, attributes>::const_iterator it = in_attributes.find(parameter);
	if (it==in_attributes.end()) return false;
	
490
	const bool isPrecip = (parameter==MeteoGrids::PSUM || parameter==MeteoGrids::SWE);
491
	const bool status = read2DGrid_internal(grid_out, filename, it->second.var, date, isPrecip);
492
493
494
495
	
	return status;
}

496
497
//isPrecip is used to convert the units of precip as well as prevent very small precip amounts
bool NetCDFIO::read2DGrid_internal(Grid2DObject& grid_out, const std::string& filename, const std::string& varname, const Date& date, const bool& isPrecip)
498
{
499
	int ncid, varid;
500
	vector<int> dimid, dim_varid;
501
502
503
	vector<string> dimname;
	vector<size_t> dimlen;

504
	if (!IOUtils::fileExists(filename)) throw AccessException(filename, AT); //prevent invalid filenames
505
	ncpp::open_file(filename, NC_NOWRITE, ncid);
506
	if (!ncpp::check_variable(ncid, varname)) return false;
507
508
	ncpp::get_variable(ncid, varname, varid);
	ncpp::get_dimension(ncid, varname, varid, dimid, dim_varid, dimname, dimlen);
509

510
	size_t time_index = IOUtils::npos, lat_index = IOUtils::npos, lon_index = IOUtils::npos;
511
	for (size_t ii=0; ii<dimname.size(); ii++) {
512
513
514
515
516
517
518
519
520
		const string name=dimname[ii];
		if (name=="latitude" || name=="lat")
			lat_index = ii;
		else if (name=="longitude" || name=="lon")
			lon_index = ii;
		else if (name=="time")
			time_index = ii;
		else
			throw InvalidArgumentException("Unknown dimension \'"+name+"\' found in file \'"+filename+"\'", AT);
521
	}
522
523
524
525
526
	if (lat_index==IOUtils::npos || lon_index==IOUtils::npos)
		throw InvalidArgumentException("Both latitudes and longitudes must be provided in file \'"+filename+"\'!", AT);
	if (dimlen[lat_index]<2 || dimlen[lon_index]<2)
		throw IOException("All dimensions for variable '" + varname + "' have to at least have length 2", AT);
	
527
	//read latitude and longitude vectors
528
529
	double *lat = new double[dimlen[lat_index]];
	double *lon = new double[dimlen[lon_index]];
530
531
	ncpp::read_data(ncid, dimname[lat_index], dim_varid[lat_index], lat);
	ncpp::read_data(ncid, dimname[lon_index], dim_varid[lon_index], lon);
532
	
533
	//read gridded data
534
535
536
	const bool date_requested = (date!=Date());
	const bool date_in_file = (time_index!=IOUtils::npos);
	
537
	double *grid = new double[dimlen[lat_index]*dimlen[lon_index]];
538
539
540
541
542
543
544
	if (!date_in_file) {
		if (!date_requested)
			ncpp::read_data(ncid, varname, varid, grid);
		else
			throw InvalidArgumentException("File \'"+filename+"\' does not contain a time dimension", AT);
	} else {
		if (date_requested) {
Mathias Bavay's avatar
Mathias Bavay committed
545
			getTimeTransform(ncid, in_time_offset, in_time_multiplier);
546
			const double timestamp = (date.getJulian() - in_time_offset) / in_time_multiplier;
547
			const size_t pos = ncpp::find_record(ncid, NetCDFIO::cf_time, timestamp);
548
549
			if (pos == IOUtils::npos) {
				double min, max;
550
				const bool status = ncpp::get_dimensionMinMax(ncid, NetCDFIO::cf_time, min, max);
551
552
553
554
				if (status) {
					Date d_min, d_max;
					d_min.setDate(min*in_time_multiplier + in_time_offset, in_dflt_TZ);
					d_max.setDate(max*in_time_multiplier + in_time_offset, in_dflt_TZ);
Mathias Bavay's avatar
Mathias Bavay committed
555
					throw IOException("No record for date " + date.toString(Date::ISO) + ". Records between " + d_min.toString(Date::ISO) + " and " + d_max.toString(Date::ISO)+" in file '"+filename+"'", AT);
556
557
558
				}
				else 
					throw IOException("No record for date " + date.toString(Date::ISO), AT);
559
			}
560
561
			ncpp::read_data(ncid, varname, varid, pos, dimlen[lat_index], dimlen[lon_index], grid);
		} else {
562
			const size_t pos = 0;
563
			ncpp::read_data(ncid, varname, varid, pos, dimlen[lat_index], dimlen[lon_index], grid);
564
		}
565
	}
566
	
567
	//read nodata value
568
	double missing_value=plugin_nodata;
569
	if (ncpp::check_attribute(ncid, varid, "missing_value")) ncpp::get_attribute(ncid, varname, varid, "missing_value", missing_value);
570

571
	//fill our Grid2DObject with all the data that has been read
572
	ncpp::copy_grid(coordin, coordinparam, dimlen[lat_index], dimlen[lon_index], lat, lon, grid, missing_value, grid_out);
573
	delete[] lat; delete[] lon; delete[] grid;
574
575
576
577
578

	//handle data packing if necessary
	if (ncpp::check_attribute(ncid, varid, "scale_factor")) {
		double scale_factor=1.;
		ncpp::get_attribute(ncid, varname, varid, "scale_factor", scale_factor);
579
		grid_out *= scale_factor;
580
581
582
583
	}
	if (ncpp::check_attribute(ncid, varid, "add_offset")) {
		double add_offset=0.;
		ncpp::get_attribute(ncid, varname, varid, "add_offset", add_offset);
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
		grid_out += add_offset;
	}
	
	//Correct the units if necessary
	if (ncpp::check_attribute(ncid, varid, "units")) {
		string units;
		ncpp::get_attribute(ncid, varname, varid, "units", units);
		if (units=="m**2 s**-2") grid_out /= Cst::gravity;
		if (units=="%") grid_out /= 100.;
		if (units=="J m**-2") grid_out /= (3600.*3.);
		if (isPrecip && units=="m") grid_out *= 100.;
		if (isPrecip) {//reset very low precip to zero
			for (size_t ii=0; ii<(grid_out.getNx()*grid_out.getNy()); ii++)
				if (grid_out(ii)<1e-3) grid_out(ii)=0.;
		}
599
	}
600

601
	ncpp::close_file(filename, ncid);
602
	return true;
603
604
}

605
void NetCDFIO::write2DGrid_internal(Grid2DObject grid_in, const std::string& filename, const attributes& attr, const Date& date, const bool& isPrecip)
606
{
607
608
	const string varname = attr.var;
	const bool is_record = (date != Date() && date!=Date(0.));
609

610
611
	double *lat_array = new double[grid_in.getNy()];
	double *lon_array = new double[grid_in.getNx()];
612
	int *data = new int[grid_in.getNy() * grid_in.getNx()];
613
614
615
	
	//Correct the units if necessary
	const string units = attr.units;
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
	if (units=="m**2 s**-2") grid_in *= Cst::gravity;
	if (units=="%") grid_in *= 100.;
	if (units=="J m**-2") grid_in *= (3600.*1.); //HACK: assuming that we do hourly outputs
	if (isPrecip && units=="m") grid_in *= 0.01;
	
	//Compute data packing
	const double data_min = grid_in.grid2D.getMin();
	const double range = grid_in.grid2D.getMax() - data_min;
	double add_offset = 0., scale_factor = 1.;
	if (range!=0.) {
		const long double type_range = 0.5* (static_cast<long double>(std::numeric_limits<int>::max()) - static_cast<long double>(std::numeric_limits<int>::min()));
		scale_factor = static_cast<double> ( range / (type_range - 1.) ); //we reserve the max for nodata
		add_offset = data_min + range/2.; //center the data on the central value of the type
		grid_in -= add_offset;
		grid_in /= scale_factor;
631
	}
632
633
634
	const double nodata_out = (range!=0)? std::numeric_limits<int>::max() : IOUtils::nodata;
	ncpp::calculate_dimensions(grid_in, lat_array, lon_array);
	ncpp::fill_grid_data(grid_in, nodata_out, data);
635
636
637
638

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

639
	if ( IOUtils::fileExists(filename) ) {
640
		ncpp::open_file(filename, NC_WRITE, ncid);
641

642
		//check of lat/lon are defined and consistent
643
		if (ncpp::check_dim_var(ncid, cf_latitude) && ncpp::check_dim_var(ncid, cf_longitude)) {
644
645
646
647
			check_consistency(ncid, grid_in, lat_array, lon_array, did_lat, did_lon, vid_lat, vid_lon);
		} else {
			create_dimensions = true;
		}
648

649
650
		if (is_record) {
			//check if a time dimension/variable already exists
651
652
653
			if (ncpp::check_dim_var(ncid, NetCDFIO::cf_time)) {
				ncpp::get_dimension(ncid, NetCDFIO::cf_time, did_time);
				ncpp::get_variable(ncid, NetCDFIO::cf_time, vid_time);
654
655
656
			} else {
				create_time = true;
			}
657
		}
658

659
660
		if (ncpp::check_variable(ncid, attr.var)) { // variable exists
			ncpp::get_variable(ncid, attr.var, vid_var);
661

662
663
664
			vector<int> dimid, dim_varid;
			vector<string> dimname;
			vector<size_t> dimlen;
665
			ncpp::get_dimension(ncid, attr.var, vid_var, dimid, dim_varid, dimname, dimlen);
666

667
			if (is_record) {
668
				if ((dimname.size() != 3) || (dimname[0] != cf_time) || (dimname[1] != cf_latitude) || (dimname[2] != cf_longitude) || (dimlen[1]!=grid_in.getNy()) || (dimlen[2]!=grid_in.getNx()))
669
					throw IOException("Variable '" + attr.var  + "' already defined with different dimensions in file '"+ filename  +"'", AT);
670
			} else {
671
				if ((dimname[0] != cf_latitude) || (dimname[1] != cf_longitude) || (dimlen[0]!=grid_in.getNy()) || (dimlen[1]!=grid_in.getNx()))
672
					throw IOException("Variable '" + attr.var  + "' already defined with different dimensions in file '"+ filename  +"'", AT);
673
			}
674
675
676
677
		} else {
			create_variable = true;
		}

678
		ncpp::start_definitions(filename, ncid);
679
	} else {
680
		if (!IOUtils::validFileAndPath(filename)) throw InvalidNameException(filename, AT);
681
682
		ncpp::create_file(filename, NC_CLASSIC_MODEL, ncid);
		ncpp::add_attribute(ncid, NC_GLOBAL, "Conventions", "CF-1.3");
683
684
		create_variable = create_dimensions = true;
		if (is_record) create_time = true;
685
686
687
688
689
	}

	if (create_dimensions) create_latlon_dimensions(ncid, grid_in, did_lat, did_lon, vid_lat, vid_lon);
	if (create_time) create_time_dimension(ncid, did_time, vid_time);

690
	if (is_record && create_variable) {
691
692
		ncpp::add_3D_variable(ncid, attr.var, NC_INT, did_time, did_lat, did_lon, vid_var); //NC_DOUBLE or NC_INT or NC_SHORT
		add_attributes_for_variable(ncid, vid_var, attr, nodata_out);
693
	} else if (create_variable) {
694
695
		ncpp::add_2D_variable(ncid, attr.var, NC_INT, did_lat, did_lon, vid_var); //NC_DOUBLE
		add_attributes_for_variable(ncid, vid_var, attr, nodata_out);
696
697
	}

698
699
700
701
702
	if (range!=0.) {
		ncpp::add_attribute(ncid, vid_var, "add_offset", add_offset);
		ncpp::add_attribute(ncid, vid_var, "scale_factor", scale_factor);
	}
	
703
	ncpp::end_definitions(filename, ncid);
704
705

	if (create_dimensions) {
706
707
		ncpp::write_data(ncid, cf_latitude, vid_lat, lat_array);
		ncpp::write_data(ncid, cf_longitude, vid_lon, lon_array);
708
709
	}

710
	if (is_record) {
711
		size_t pos_start = ncpp::add_record(ncid, NetCDFIO::cf_time, vid_time, static_cast<double>( date.getUnixDate()/3600 ));
712
		ncpp::write_data(ncid, attr.var, vid_var, grid_in.getNy(), grid_in.getNx(), pos_start, data);
713
	} else {
714
		ncpp::write_data(ncid, attr.var, vid_var, data);
715
	}
716

717
	ncpp::close_file(filename, ncid);
718
719
720
	delete[] lat_array; delete[] lon_array; delete[] data;
}

721
void NetCDFIO::getTimeTransform(const int& ncid, double &time_offset, double &time_multiplier) const
722
723
{
	string time_units;
724
725
	ncpp::get_DimAttribute(ncid, NetCDFIO::cf_time, "units", time_units);

726
727
728
729
	std::vector<std::string> vecString;
	const size_t nrWords = IOUtils::readLineToVec(time_units, vecString);
	if (nrWords<3 || nrWords>4) throw InvalidArgumentException("Invalid format for time units: \'"+time_units+"\'", AT);
	
730
731
732
733
734
	const double equinox_year = 365.242198781; //definition used by the NetCDF Udunits package
	
	if (vecString[0]=="years") time_multiplier = equinox_year;
	else if (vecString[0]=="months") time_multiplier = equinox_year/12.;
	else if (vecString[0]=="days") time_multiplier = 1.;
735
736
737
738
739
740
741
742
743
744
745
746
747
	else if (vecString[0]=="hours") time_multiplier = 1./24.;
	else if (vecString[0]=="minutes") time_multiplier = 1./(24.*60.);
	else if (vecString[0]=="seconds") time_multiplier = 1./(24.*3600);
	else throw InvalidArgumentException("Unknown time unit \'"+vecString[0]+"\'", AT);
	
	const string ref_date_str = (nrWords==3)? vecString[2] : vecString[2]+"T"+vecString[3];
	Date refDate;
	if (!IOUtils::convertString(refDate, ref_date_str, in_dflt_TZ))
		throw InvalidArgumentException("Invalid reference date \'"+ref_date_str+"\'", AT);
	
	time_offset = refDate.getJulian();
}

748
void NetCDFIO::create_latlon_dimensions(const int& ncid, const Grid2DObject& grid_in, int& did_lat, int& did_lon, int& varid_lat, int& varid_lon)
749
{
750
	ncpp::add_dimension(ncid, cf_latitude, grid_in.getNy(), did_lat);
751
752
753
754
	ncpp::add_1D_variable(ncid, cf_latitude, NC_DOUBLE, did_lat, varid_lat);
	ncpp::add_attribute(ncid, varid_lat, "standard_name", cf_latitude);
	ncpp::add_attribute(ncid, varid_lat, "long_name", cf_latitude);
	ncpp::add_attribute(ncid, varid_lat, "units", "degrees_north");
755

756
	ncpp::add_dimension(ncid, cf_longitude, grid_in.getNx(), did_lon);
757
758
759
760
	ncpp::add_1D_variable(ncid, cf_longitude, NC_DOUBLE, did_lon, varid_lon);
	ncpp::add_attribute(ncid, varid_lon, "standard_name", cf_longitude);
	ncpp::add_attribute(ncid, varid_lon, "long_name", cf_longitude);
	ncpp::add_attribute(ncid, varid_lon, "units", "degrees_east");
761
762
}

763
void NetCDFIO::create_time_dimension(const int& ncid, int& did_time, int& varid_time)
764
{
765
	ncpp::add_dimension(ncid, NetCDFIO::cf_time, NC_UNLIMITED, did_time);
766
767
768
	ncpp::add_1D_variable(ncid, NetCDFIO::cf_time, NC_DOUBLE, did_time, varid_time);
	ncpp::add_attribute(ncid, varid_time, "standard_name", cf_time);
	ncpp::add_attribute(ncid, varid_time, "long_name", cf_time);
769
	ncpp::add_attribute(ncid, varid_time, "units", "hours since 1970-1-1");
770
	ncpp::add_attribute(ncid, varid_time, "calendar", "gregorian");
771
772
}

773
void NetCDFIO::add_attributes_for_variable(const int& ncid, const int& varid, const attributes& attr, const double& nodata_out)
774
{
775
776
777
778
779
780
	ncpp::add_attribute(ncid, varid, "standard_name", attr.standard_name);
	ncpp::add_attribute(ncid, varid, "long_name", attr.long_name);
	ncpp::add_attribute(ncid, varid, "units", attr.units);
	if (attr.var=="z" || attr.var=="ZS") { //for DEM
		ncpp::add_attribute(ncid, varid, "positive", "up");
		ncpp::add_attribute(ncid, varid, "axis", "Z");
781
	}
782
	ncpp::add_attribute(ncid, varid, "missing_value", nodata_out);
783
}
784
785
786
787
788
789

void NetCDFIO::check_consistency(const int& ncid, const Grid2DObject& grid, double*& lat_array, double*& lon_array,
                                 int& did_lat, int& did_lon, int& vid_lat, int& vid_lon)
{
	size_t latlen, lonlen;

790
791
	ncpp::get_dimension(ncid, cf_latitude, did_lat, latlen);
	ncpp::get_dimension(ncid, cf_longitude, did_lon, lonlen);
792

793
794
	ncpp::get_variable(ncid, cf_latitude, vid_lat);
	ncpp::get_variable(ncid, cf_longitude, vid_lon);
795

796
	if ((latlen != grid.getNy()) || (lonlen != grid.getNx()))
797
798
		throw IOException("Error while writing grid - grid size and lat/lon coordinates are inconsistent", AT);

799
800
	double *lat = new double[grid.getNy()];
	double *lon = new double[grid.getNx()];
801

802
803
	ncpp::read_data(ncid, cf_latitude, vid_lat, lat);
	ncpp::read_data(ncid, cf_longitude, vid_lon, lon);
804

805
	for (size_t ii=0; ii<latlen; ++ii) {
806
807
808
809
		if (lat_array[ii] != lat[ii])
			throw IOException("Error while writing grid - grid and lat/lon coordinates are inconsistent", AT);
	}

810
	for (size_t ii=0; ii<lonlen; ++ii) {
811
812
813
		if (lon_array[ii] != lon[ii])
			throw IOException("Error while writing grid - grid and lat/lon coordinates are inconsistent", AT);
	}
Thomas Egger's avatar
Thomas Egger committed
814
815

	delete[] lat; delete[] lon;
816
817
}

818
} //namespace