WSL/SLF GitLab Repository

NetCDFIO.cc 39.9 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
44
45
46
 * <A HREF="http://www.epic.noaa.gov/java/ncBrowse/">ncBrowse</A> java software or 
 * <A HREF="http://meteora.ucsd.edu/~pierce/ncview_home_page.html">ncview</A>. It is also possible to run *ncdump* on a given
 * file in order to have a look at its structure (such as *ncdump {my_netcdf_file} | more*) and specially the parameters names 
 * (this is useful if remapping is needed, see below for in the \ref netcdf_keywords "keywords" section).
47
 *
48
49
50
51
 * 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
52
 * - 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;
53
 * - CNRM - from the <A HREF="http://www.cnrm.meteo.fr/">National Centre for Meteorological Research</A>.
54
55
56
57
 * 
 * @section netcdf_compilation Compilation
 * 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.
58
59
60
61
62
 *
 * @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
63
64
 * - 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
65
66
67
68
 * - 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
69
 * - NETCDF_SCHEMA: the schema to use (either CF1 or CNRM or ECMWF); [Input] and [Output] section
70
 * - 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;
71
 * - 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
72
 *
Mathias Bavay's avatar
Mathias Bavay committed
73
74
75
 * 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.
 * 
76
 * @section netcdf_example Example use
77
 * Using this plugin to build downscaled time series at virtual stations, with the ECMWF Era Interim data set (see section below):
78
79
 * @code
 * [Input]
Mathias Bavay's avatar
Mathias Bavay committed
80
81
82
83
 * GRID2D    = NETCDF
 * GRID2DPATH =  /data/meteo_reanalysis
 * METEO_EXT = .nc
 * NETCDF_SCHEMA = ECMWF
84
85
86
87
88
89
90
91
92
 * 
 * DEM = NETCDF
 * DEMFILE = /data/meteo_reanalysis/ECMWF_Europe_20150101-20150701.nc
 * DEM_FROM_PRESSURE = true
 * 
 * #The lines below have nothing to do with this plugin
 * Downscaling = true
 * VSTATION1 = 46.793029 9.821343 ;this is Davos
 * Virtual_parameters = TA RH PSUM ISWR ILWR P VW DW TSS HS RSWR TSG ;this has to fit the parameter set in the data files
93
 * @endcode
94
 * 
95
 * Another example, to extract precipitation from the MeteoSwiss daily precipitation reanalysis, RhiresD
96
 * @code
97
98
99
100
 * [Input]
 * DEM     = NETCDF
 * DEMFILE = ./input/ch02_lonlat.nc
 * 
101
102
103
 * GRID2D    = NETCDF
 * GRID2DPATH =  /data/meteo_reanalysis
 * METEO_EXT = .nc
104
 * NETCDF::PSUM = RhiresD               ;overwrite the PSUM parameter with "RhiresD", for example for MeteoCH reanalysis
105
 * 
106
 * #The lines below have nothing to do with this plugin
107
108
 * Downscaling = true
 * VSTATION1 = 46.793029 9.821343 ;this is Davos
109
 * Virtual_parameters = PSUM ;this has to fit the parameter set in the data files
110
 * @endcode
111
112
113
 * 
 * @section netcdf_meteoch MeteoCH RhiresD & similar products
 * <A HREF="http://www.meteoswiss.admin.ch/home.html?tab=overview">MeteoSwiss</A> provides <A HREF="http://www.ifu.ethz.ch/hydrologie/research/research_data/proddoc.pdf">reanalysis</A> of precipitation and other meteo fields from 1961 to present over Switzerland for different time scales: daily, monthly, yearly, as well as hourly (CPC dataset). The DEM are also provided, either in lat/lon, 
114
 * Swiss coordinates, rotated lat/lon, ... These data sets must be requested from MeteoSwiss and are available with a specific license for research. 
115
 *
116
117
118
119
120
121
122
123
124
 * Unfortunatelly, the naming of the parameters and dimensions within the files is not standard nor consistent. In order to handle the parameters names, 
 * simply run *ncdump {my_netcdf_file} | more* and use the name mapping facility of this plugin to map the non-standard parameters to our internal names
 * (see the \ref netcdf_keywords "plugin keywords"). When the dimensions are not standard (for example the time axis being called "TS_REFERENCE"), 
 * use first the <A HREF="http://linux.die.net/man/1/ncrename">ncrename</A> tool that is part of the 
 * <A HREF="http://nco.sourceforge.net/">NCO utilities</A> to rename both the dimension (-d) and the variable (-v): 
 * @code
 * ncrename -d TS_REFERENCE,time -v TS_REFERENCE,time {my_netcdf_file}
 * @endcode
 * 
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
 * @section netcdf_ecmwf ECMWF Era Interim
 * The Era Interim data can be downloaded on the <A HREF="http://apps.ecmwf.int/datasets/data/interim-full-daily/levtype=sfc/">ECMWF dataserver</A> 
 * after creating an account and login in. 
 * 
 * It is recommended to extract data at 00:00, and 12:00 for all steps 3, 6, 9, 12. The select the following fields:
 * 10 metre U wind component, 10 metre V wind component, 2 metre dewpoint temperature, 2 metre temperature, Forecast albedo, Mean sea level pressure, Skin temperature, Snow density, Snow depth, Soil temperature level 1, Surface pressure, Surface solar radiation downwards, Surface thermal radiation downwards, Total precipitation
 * 
 * Here we have included the *forecast albedo* so the RSWR can be computed from ISWR and the *mean sea level pressure* and *surface pressure*
 * as proxies to compute the elevation. If you have the altitude in a separate file, it can be declared as DEM and there would be no need for the sea 
 *level pressure (this would also be much more precise).
 * 
 * You should therefore have the following request:
 * @code
 * Parameter: 10 metre U wind component, 10 metre V wind component, 2 metre dewpoint temperature, 2 metre temperature, Forecast albedo, 
 *            Mean sea level pressure, Skin temperature, Snow density, Snow depth, Soil temperature level 1, Surface pressure, 
 *            Surface solar radiation downwards, Surface thermal radiation downwards, Total precipitation
 *      Step: 3 to 12 by 3
 *      Type: Forecast
 *      Time: 00:00:00, 12:00:00
 * @endcode
 * 
 * With the <A HREF="https://software.ecmwf.int/wiki/display/WEBAPI/Access+ECMWF+Public+Datasets">ECMWF Python Library</A>, the request 
 * would be for example:
 * @code
 * #!/usr/bin/env python
 * from ecmwfapi import ECMWFDataServer
 * server = ECMWFDataServer()
 * server.retrieve({
153
154
155
156
157
158
159
160
161
162
163
164
165
166
 * "class": "ei",
 * "dataset": "interim",
 * "date": "2015-01-01/to/2015-01-31",
 * "expver": "1",
 * "grid": "0.75/0.75",
 * "levtype": "sfc",
 * "param": "33.128/134.128/139.128/141.128/151.128/165.128/166.128/167.128/168.128/169.128/175.128/205.128/228.128/235.128/243.128",
 * "step": "3/6/9/12",
 * "area":"42.2/-1.5/51.7/15.7",
 * "stream": "oper",
 * "format":"netcdf",
 * "target": "my-era-interim.nc",
 * "time": "00/12",
 * "type": "fc",
167
 * })
168
 * @endcode
169
 */
170
 
171
const double NetCDFIO::plugin_nodata = -9999999.; //CNRM-GAME nodata value
172
const std::string NetCDFIO::cf_time = "time";
173
174
const std::string NetCDFIO::cf_latitude = "latitude";
const std::string NetCDFIO::cf_longitude = "longitude";
175

176
177
NetCDFIO::NetCDFIO(const std::string& configfile) : cfg(configfile), cache_meteo_files(), in_attributes(), out_attributes(), 
                                                    coordin(), coordinparam(), coordout(), coordoutparam(), 
178
                                                    in_dflt_TZ(0.), out_dflt_TZ(0.), in_time_offset(IOUtils::nodata), in_time_multiplier(IOUtils::nodata), 
179
                                                    dem_altimeter(false), in_strict(false), out_strict(false), vecMetaData()
180
181
{
	IOUtils::getProjectionParameters(cfg, coordin, coordinparam, coordout, coordoutparam);
182
	parseInputOutputSection();
183
184
}

185
186
NetCDFIO::NetCDFIO(const Config& cfgreader) : cfg(cfgreader), cache_meteo_files(), in_attributes(), out_attributes(), 
                                              coordin(), coordinparam(), coordout(), coordoutparam(), 
187
                                              in_dflt_TZ(0.), out_dflt_TZ(0.), in_time_offset(IOUtils::nodata), in_time_multiplier(IOUtils::nodata), 
188
                                              dem_altimeter(false), in_strict(false), out_strict(false), vecMetaData()
189
190
{
	IOUtils::getProjectionParameters(cfg, coordin, coordinparam, coordout, coordoutparam);
191
	parseInputOutputSection();
192
193
}

194
195
196
197
198
199
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);
200
	cfg.getValue("DEM_FROM_PRESSURE", "Input", dem_altimeter, IOUtils::nothrow);
201
	
202
203
204
205
206
207
208
	const string in_schema = IOUtils::strToUpper( cfg.get("NETCDF_SCHEMA", "Input", IOUtils::nothrow) );
	if (!in_schema.empty()) initAttributesMap(in_schema, in_attributes);
	else initAttributesMap("ECMWF", in_attributes);
	
	const string out_schema = IOUtils::strToUpper( cfg.get("NETCDF_SCHEMA", "Output", IOUtils::nothrow) );
	if (!out_schema.empty()) initAttributesMap(out_schema, out_attributes);
	else initAttributesMap("ECMWF", out_attributes);
209
	
210
	const std::string in_meteo = cfg.get("GRID2D", "Input", IOUtils::nothrow);
211
212
213
214
	if (in_meteo == "NETCDF") { //keep it synchronized with IOHandler.cc for plugin mapping!!
		const string in_grid2d_path = cfg.get("GRID2DPATH", "Input", IOUtils::nothrow);
		if (!in_grid2d_path.empty()) scanMeteoPath(in_grid2d_path,  cache_meteo_files);
	}
215
216
}

217
void NetCDFIO::initAttributesMap(const std::string& schema, std::map<MeteoGrids::Parameters, attributes> &attr)
218
219
220
221
222
223
224
225
226
227
228
229
230
231
{
	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);
232
		attr[MeteoGrids::QI] = attributes("Qair", "", "", "", IOUtils::nodata);
233
234
		attr[MeteoGrids::VW] = attributes("Wind", "", "Wind Speed", "m/s", IOUtils::nodata);
		attr[MeteoGrids::DW] = attributes("Wind_DIR", "", "Wind Direction", "deg", IOUtils::nodata);
235
236
		attr[MeteoGrids::PSUM_L] = attributes("Rainf", "", "Rainfall Rate", "kg/m2/s", IOUtils::nodata);
		attr[MeteoGrids::PSUM_S] = attributes("Snowf", "", "", "", IOUtils::nodata);
237
238
		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);
239
240
241
		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") {
242
		attr[MeteoGrids::DEM] = attributes("z", "geopotential_height", "geopotential_height", "m", IOUtils::nodata);
243
244
245
246
247
248
		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);
249
		attr[MeteoGrids::PSUM] = attributes("tp", "", "Total precipitation", "m", IOUtils::nodata);
250
251
252
253
		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);
254
		attr[MeteoGrids::TSG] = attributes("stl1", "surface_temperature", "Soil temperature level 1", "K", IOUtils::nodata); //this is from 0 to -7cm
255
256
		//attr[MeteoGrids::ALB] = attributes("al", "surface_albedo", "Albedo", "(0 - 1)", IOUtils::nodata);
		attr[MeteoGrids::ALB] = attributes("fal", "", "Forecast albedo", "(0 - 1)", IOUtils::nodata);
257
258
		attr[MeteoGrids::RSNO] = attributes("rsn", "", "Snow density", "kg m**-3", IOUtils::nodata);
		attr[MeteoGrids::ROT] = attributes("ro", "", "Runoff", "m", IOUtils::nodata);
259
260
	} else
		throw InvalidArgumentException("Invalid schema selected for NetCDF: \""+schema+"\"", AT);
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
	
	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);
	}
276
277
}

278
void NetCDFIO::read2DGrid(Grid2DObject& grid_out, const std::string& arguments)
279
{
280
281
282
283
284
285
286
287
	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);
	}
288
289
}

290
void NetCDFIO::read2DGrid(Grid2DObject& grid_out, const MeteoGrids::Parameters& parameter, const Date& date)
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
{
	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)
313
{
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
	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) {
345
			for (size_t ii=0; ii<(grid_out.getNx()*grid_out.getNy()); ii++)
346
347
348
349
350
351
				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) {
352
			for (size_t ii=0; ii<(grid_out.getNx()*grid_out.getNy()); ii++)
353
354
355
356
357
				grid_out(ii) = Atmosphere::DewPointtoRh(grid_out(ii), ta(ii), false);
			return;
		}
	}
	
358
359
360
361
362
363
364
365
366
367
368
369
370
371
	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;
		}
	}
	
372
373
374
375
376
377
378
379
380
381
	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;
		}
	}
	
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
	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;
			}
404
405
406
407
			return;
		}
	}
	
408
409
410
411
412
413
414
415
416
417
418
	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;
		}
	}
	
419
	throw InvalidArgumentException("Parameter \'"+MeteoGrids::getParameterName(parameter)+"\' either not found in file \'"+filename+"\' or not found in current NetCDF schema", AT);
420
421
}

422
void NetCDFIO::readDEM(DEMObject& dem_out)
423
{
424
425
426
427
428
429
	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 {
430
		if (read2DGrid_internal(dem_out, filename, MeteoGrids::DEM)) return; //schema naming
431
432
		if (read2DGrid_internal(dem_out, filename, "Band1")) return; //ASTER naming
		if (read2DGrid_internal(dem_out, filename, "z")) return; //GDAL naming
433
		if (read2DGrid_internal(dem_out, filename, "height")) return; //MeteoCH naming
434
		
435
		//last chance: read from pressure grids
436
437
		if (dem_altimeter) {
			Grid2DObject p, ta, p_sea;
438
439
440
			if (read2DGrid_internal(p_sea, filename, MeteoGrids::P_SEA) &&
			read2DGrid_internal(p, filename, MeteoGrids::P) && 
			read2DGrid_internal(ta, filename, MeteoGrids::TA)) {
441
				dem_out.set(p, IOUtils::nodata);
442
				const double k = Cst::gravity / (Cst::mean_adiabatique_lapse_rate * Cst::gaz_constant_dry_air);
443
				const double k_inv = 1./k;
444
				for (size_t ii=0; ii<(dem_out.getNx()*dem_out.getNy()); ii++) {
445
					const double K = pow(p(ii)/p_sea(ii), k_inv);
446
					dem_out(ii) = ta(ii)*Cst::earth_R0*(1.-K) / (Cst::mean_adiabatique_lapse_rate * Cst::earth_R0 - ta(ii)*(1.-K));
447
				}
448
449
				return;
			}
450
451
		}
		
452
453
454
455
456
457
458
459
		throw InvalidArgumentException("The variable containing the DEM could not be found. Please specify it using the DEMVAR key.", 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;
460
	if (IOUtils::readLineToVec(arguments, vec_argument, ':')  != 2)
461
		throw InvalidArgumentException("The format for the arguments to NetCDFIO::write2DGrid is filename:varname", AT);
462

463
464
465
	const string name = vec_argument[1];
	const attributes attr(name, name, name, "", IOUtils::nodata);
	write2DGrid_internal(grid_in, vec_argument[0], attr);
466
467
468
469
470
}

void NetCDFIO::write2DGrid(const Grid2DObject& grid_in, const MeteoGrids::Parameters& parameter, const Date& date)
{
	const string filename = cfg.get("GRID2DFILE", "Output");
471
472
473
	
	const std::map<MeteoGrids::Parameters, attributes>::const_iterator it = in_attributes.find(parameter);
	if (it!=in_attributes.end()) {
474
		const bool isPrecip = (parameter==MeteoGrids::PSUM || parameter==MeteoGrids::SWE);
475
476
477
478
479
480
		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);
	}
481
482
}

483
484
485
486
487
488
489
490
491
492
493
494
495
//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();

496
	const string meteo_ext = cfg.get("METEO_EXT", "INPUT", IOUtils::nothrow);
497
	std::list<std::string> dirlist = FileUtils::readDirectory(meteopath_in, meteo_ext);
498
499
500
501
502
503
504
505
	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;
		
506
		if (!FileUtils::fileExists(filename)) throw AccessException(filename, AT); //prevent invalid filenames
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
		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;
	}
}

532
533
534
535
536
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;
	
537
	const bool isPrecip = (parameter==MeteoGrids::PSUM || parameter==MeteoGrids::SWE);
538
	const bool status = read2DGrid_internal(grid_out, filename, it->second.var, date, isPrecip);
539
540
541
542
	
	return status;
}

543
544
//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)
545
{
546
	int ncid, varid;
547
	vector<int> dimid, dim_varid;
548
549
550
	vector<string> dimname;
	vector<size_t> dimlen;

551
	if (!FileUtils::fileExists(filename)) throw AccessException(filename, AT); //prevent invalid filenames
552
	ncpp::open_file(filename, NC_NOWRITE, ncid);
553
	if (!ncpp::check_variable(ncid, varname)) return false;
554
555
	ncpp::get_variable(ncid, varname, varid);
	ncpp::get_dimension(ncid, varname, varid, dimid, dim_varid, dimname, dimlen);
556

557
	size_t time_index = IOUtils::npos, lat_index = IOUtils::npos, lon_index = IOUtils::npos;
558
	for (size_t ii=0; ii<dimname.size(); ii++) {
559
560
561
562
563
564
565
566
567
		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);
568
	}
569
570
571
572
573
	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);
	
574
	//read latitude and longitude vectors
575
576
	double *lat = new double[dimlen[lat_index]];
	double *lon = new double[dimlen[lon_index]];
577
578
	ncpp::read_data(ncid, dimname[lat_index], dim_varid[lat_index], lat);
	ncpp::read_data(ncid, dimname[lon_index], dim_varid[lon_index], lon);
579
	
580
	//read gridded data
581
582
583
	const bool date_requested = (date!=Date());
	const bool date_in_file = (time_index!=IOUtils::npos);
	
584
	double *grid = new double[dimlen[lat_index]*dimlen[lon_index]];
585
586
587
588
589
590
591
	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
592
			getTimeTransform(ncid, in_time_offset, in_time_multiplier);
593
			const double timestamp = (date.getJulian() - in_time_offset) / in_time_multiplier;
594
			const size_t pos = ncpp::find_record(ncid, NetCDFIO::cf_time, timestamp);
595
596
			if (pos == IOUtils::npos) {
				double min, max;
597
				const bool status = ncpp::get_dimensionMinMax(ncid, NetCDFIO::cf_time, min, max);
598
599
600
601
				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
602
					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);
603
604
605
				}
				else 
					throw IOException("No record for date " + date.toString(Date::ISO), AT);
606
			}
607
608
			ncpp::read_data(ncid, varname, varid, pos, dimlen[lat_index], dimlen[lon_index], grid);
		} else {
609
			const size_t pos = 0;
610
			ncpp::read_data(ncid, varname, varid, pos, dimlen[lat_index], dimlen[lon_index], grid);
611
		}
612
	}
613
	
614
	//read nodata value
615
	double missing_value=plugin_nodata;
616
	if (ncpp::check_attribute(ncid, varid, "missing_value")) ncpp::get_attribute(ncid, varname, varid, "missing_value", missing_value);
617

618
	//fill our Grid2DObject with all the data that has been read
619
	ncpp::copy_grid(coordin, coordinparam, dimlen[lat_index], dimlen[lon_index], lat, lon, grid, missing_value, grid_out);
620
	delete[] lat; delete[] lon; delete[] grid;
621
622
623
624
625

	//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);
626
		grid_out *= scale_factor;
627
628
629
630
	}
	if (ncpp::check_attribute(ncid, varid, "add_offset")) {
		double add_offset=0.;
		ncpp::get_attribute(ncid, varname, varid, "add_offset", add_offset);
631
632
633
634
635
636
637
638
639
640
641
642
643
		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++)
644
				if (grid_out(ii)<1e-3 && grid_out(ii)!=mio::IOUtils::nodata) grid_out(ii)=0.;
645
		}
646
	}
647

648
	ncpp::close_file(filename, ncid);
649
	return true;
650
651
}

652
void NetCDFIO::write2DGrid_internal(Grid2DObject grid_in, const std::string& filename, const attributes& attr, const Date& date, const bool& isPrecip)
653
{
654
655
	const string varname = attr.var;
	const bool is_record = (date != Date() && date!=Date(0.));
656

657
658
	double *lat_array = new double[grid_in.getNy()];
	double *lon_array = new double[grid_in.getNx()];
659
	int *data = new int[grid_in.getNy() * grid_in.getNx()];
660
661
662
	
	//Correct the units if necessary
	const string units = attr.units;
663
664
665
666
667
668
	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;
	
	ncpp::calculate_dimensions(grid_in, lat_array, lon_array);
669
	ncpp::fill_grid_data(grid_in, IOUtils::nodata, data);
670
671
672
673

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

674
	if ( FileUtils::fileExists(filename) ) {
675
		ncpp::open_file(filename, NC_WRITE, ncid);
676

677
		//check of lat/lon are defined and consistent
678
		if (ncpp::check_dim_var(ncid, cf_latitude) && ncpp::check_dim_var(ncid, cf_longitude)) {
679
680
681
682
			check_consistency(ncid, grid_in, lat_array, lon_array, did_lat, did_lon, vid_lat, vid_lon);
		} else {
			create_dimensions = true;
		}
683

684
685
		if (is_record) {
			//check if a time dimension/variable already exists
686
687
688
			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);
689
690
691
			} else {
				create_time = true;
			}
692
		}
693

694
695
		if (ncpp::check_variable(ncid, attr.var)) { // variable exists
			ncpp::get_variable(ncid, attr.var, vid_var);
696

697
698
699
			vector<int> dimid, dim_varid;
			vector<string> dimname;
			vector<size_t> dimlen;
700
			ncpp::get_dimension(ncid, attr.var, vid_var, dimid, dim_varid, dimname, dimlen);
701

702
			if (is_record) {
703
				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()))
704
					throw IOException("Variable '" + attr.var  + "' already defined with different dimensions in file '"+ filename  +"'", AT);
705
			} else {
706
				if ((dimname[0] != cf_latitude) || (dimname[1] != cf_longitude) || (dimlen[0]!=grid_in.getNy()) || (dimlen[1]!=grid_in.getNx()))
707
					throw IOException("Variable '" + attr.var  + "' already defined with different dimensions in file '"+ filename  +"'", AT);
708
			}
709
710
711
712
		} else {
			create_variable = true;
		}

713
		ncpp::start_definitions(filename, ncid);
714
	} else {
715
		if (!FileUtils::validFileAndPath(filename)) throw InvalidNameException(filename, AT);
716
717
		ncpp::create_file(filename, NC_CLASSIC_MODEL, ncid);
		ncpp::add_attribute(ncid, NC_GLOBAL, "Conventions", "CF-1.3");
718
719
		create_variable = create_dimensions = true;
		if (is_record) create_time = true;
720
721
722
723
724
	}

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

725
	if (is_record && create_variable) {
726
		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
727
		add_attributes_for_variable(ncid, vid_var, attr, IOUtils::nodata);
728
	} else if (create_variable) {
729
		ncpp::add_2D_variable(ncid, attr.var, NC_INT, did_lat, did_lon, vid_var); //NC_DOUBLE
730
		add_attributes_for_variable(ncid, vid_var, attr, IOUtils::nodata);
731
	}
732
	ncpp::end_definitions(filename, ncid);
733
734

	if (create_dimensions) {
735
736
		ncpp::write_data(ncid, cf_latitude, vid_lat, lat_array);
		ncpp::write_data(ncid, cf_longitude, vid_lon, lon_array);
737
738
	}

739
	if (is_record) {
740
		size_t pos_start = ncpp::add_record(ncid, NetCDFIO::cf_time, vid_time, static_cast<double>( date.getUnixDate()/3600 ));
741
		ncpp::write_data(ncid, attr.var, vid_var, grid_in.getNy(), grid_in.getNx(), pos_start, data);
742
	} else {
743
		ncpp::write_data(ncid, attr.var, vid_var, data);
744
	}
745

746
	ncpp::close_file(filename, ncid);
747
748
749
	delete[] lat_array; delete[] lon_array; delete[] data;
}

750
void NetCDFIO::getTimeTransform(const int& ncid, double &time_offset, double &time_multiplier) const
751
752
{
	string time_units;
753
754
	ncpp::get_DimAttribute(ncid, NetCDFIO::cf_time, "units", time_units);

755
756
757
758
	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);
	
759
760
761
762
763
	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.;
764
765
766
767
768
769
770
771
772
773
774
775
776
	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();
}

777
void NetCDFIO::create_latlon_dimensions(const int& ncid, const Grid2DObject& grid_in, int& did_lat, int& did_lon, int& varid_lat, int& varid_lon)
778
{
779
	ncpp::add_dimension(ncid, cf_latitude, grid_in.getNy(), did_lat);
780
781
782
783
	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");
784

785
	ncpp::add_dimension(ncid, cf_longitude, grid_in.getNx(), did_lon);
786
787
788
789
	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");
790
791
}

792
void NetCDFIO::create_time_dimension(const int& ncid, int& did_time, int& varid_time)
793
{
794
	ncpp::add_dimension(ncid, NetCDFIO::cf_time, NC_UNLIMITED, did_time);
795
796
797
	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);
798
	ncpp::add_attribute(ncid, varid_time, "units", "hours since 1970-1-1");
799
	ncpp::add_attribute(ncid, varid_time, "calendar", "gregorian");
800
801
}

802
void NetCDFIO::add_attributes_for_variable(const int& ncid, const int& varid, const attributes& attr, const double& nodata_out)
803
{
804
805
806
807
808
809
	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");
810
	}
811
	ncpp::add_attribute(ncid, varid, "missing_value", nodata_out);
812
}
813
814
815
816
817
818

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;

819
820
	ncpp::get_dimension(ncid, cf_latitude, did_lat, latlen);
	ncpp::get_dimension(ncid, cf_longitude, did_lon, lonlen);
821

822
823
	ncpp::get_variable(ncid, cf_latitude, vid_lat);
	ncpp::get_variable(ncid, cf_longitude, vid_lon);
824

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

828
829
	double *lat = new double[grid.getNy()];
	double *lon = new double[grid.getNx()];
830

831
832
	ncpp::read_data(ncid, cf_latitude, vid_lat, lat);
	ncpp::read_data(ncid, cf_longitude, vid_lon, lon);
833

834
	for (size_t ii=0; ii<latlen; ++ii) {
835
836
837
838
		if (lat_array[ii] != lat[ii])
			throw IOException("Error while writing grid - grid and lat/lon coordinates are inconsistent", AT);
	}

839
	for (size_t ii=0; ii<lonlen; ++ii) {
840
841
842
		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
843
844

	delete[] lat; delete[] lon;
845
846
}

847
} //namespace