WSL/SLF GitLab Repository

ARPSIO.cc 16.6 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
/***********************************************************************************/
/*  Copyright 2009 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/>.
*/
18
#include <meteoio/meteolaws/Meteoconst.h> //for PI
19
20
#include <meteoio/MathOptim.h>

21
22
23
#include <string.h>
#include <algorithm>

24
25
26
27
28
29
30
31
#include "ARPSIO.h"

using namespace std;

namespace mio {
/**
 * @page arps ARPSIO
 * @section arps_format Format
32
 * This is for reading grid data in the ARPS grid format (it transparently supports both true ARPS ascii grids and grids modified by the ARPSGRID utility). DEM reading works well while reading meteo parameters might be a rough ride (since ARPS files do not always contain a consistent set of meteo fields).
33
34
 *
 * @section arps_units Units
35
 * All units are assumed to be MKSA.
36
37
38
39
40
41
 *
 * @section arps_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
 * - DEMFILE: path and file containing the DEM; [Input] section
42
43
44
45
 * - ARPS_XCOORD: x coordinate of the lower left corner of the grids; [Input] section
 * - ARPS_YCOORD: y coordinate of the lower left corner of the grids; [Input] section
 * - GRID2DPATH: path to the input directory where to find the arps files to be read as grids; [Input] section
 * - GRID2DEXT: arps file extension, or <i>none</i> for no file extension (default: .asc)
46
47
48
 */

const double ARPSIO::plugin_nodata = -999.; //plugin specific nodata value
49
const std::string ARPSIO::default_ext=".asc"; //filename extension
50

51
ARPSIO::ARPSIO(const std::string& configfile)
52
        : cfg(configfile),
53
          fin(NULL), filename(),  coordin(), coordinparam(), coordout(), coordoutparam(),
54
55
          grid2dpath_in(), ext(default_ext), dimx(0), dimy(0), dimz(0), cellsize(0.),
          xcoord(IOUtils::nodata), ycoord(IOUtils::nodata), zcoord(), is_true_arps(true)
56
57
{
	IOUtils::getProjectionParameters(cfg, coordin, coordinparam, coordout, coordoutparam);
58
	setOptions();
59
60
}

61
ARPSIO::ARPSIO(const Config& cfgreader)
62
        : cfg(cfgreader),
63
          fin(NULL), filename(),  coordin(), coordinparam(), coordout(), coordoutparam(),
64
65
          grid2dpath_in(), ext(default_ext), dimx(0), dimy(0), dimz(0), cellsize(0.),
          xcoord(IOUtils::nodata), ycoord(IOUtils::nodata), zcoord(), is_true_arps(true)
66
67
{
	IOUtils::getProjectionParameters(cfg, coordin, coordinparam, coordout, coordoutparam);
68
69
70
	setOptions();
}

71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
ARPSIO& ARPSIO::operator=(const ARPSIO& source) {
	if(this != &source) {
		fin = NULL;
		filename = source.filename;
		coordin = source.coordin;
		coordinparam = source.coordinparam;
		coordout = source.coordout;
		coordoutparam = source.coordoutparam;
		grid2dpath_in = source. grid2dpath_in;
		ext = source.ext;
		dimx = source.dimx;
		dimy = source.dimy;
		dimz = source.dimz;
		cellsize = source.cellsize;
		xcoord = source.xcoord;
		ycoord = source.ycoord;
		zcoord = source.zcoord;
		is_true_arps = source.is_true_arps;
	}
	return *this;
}

93
94
void ARPSIO::setOptions()
{
95
	string tmp;
96
	cfg.getValue("GRID2D", "Input", tmp, IOUtils::nothrow);
97
98
99
100
	if (tmp == "ARPS") { //keep it synchronized with IOHandler.cc for plugin mapping!!
		cfg.getValue("GRID2DPATH", "Input", grid2dpath_in);
	}

101
102
	cfg.getValue("ARPS_XCOORD", "Input", xcoord, IOUtils::dothrow);
	cfg.getValue("ARPS_YCOORD", "Input", ycoord, IOUtils::dothrow);
103

104
	//default value has been set in constructor
105
	cfg.getValue("GRID2DEXT", "Input", ext, IOUtils::nothrow);
106
	if(ext=="none") ext.clear();
107
108
109
110
111
112
113
}

ARPSIO::~ARPSIO() throw()
{
	cleanup();
}

114
void ARPSIO::read2DGrid(Grid2DObject& grid_out, const std::string& i_name)
115
{
116
	const std::string _filename = grid2dpath_in +"/" + i_name;
117
118

	openGridFile(_filename);
119
120
121
122
123
124

	const unsigned int layer=2;
	if(is_true_arps)
		readGridLayer("zp coordinat", layer, grid_out);
	else
		readGridLayer("zp_coordinat", layer, grid_out);
125
126
127
128
129

	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
void ARPSIO::read2DGrid(Grid2DObject& grid_out, const MeteoGrids::Parameters& parameter, const Date& date)
{
	std::string date_str = date.toString(Date::ISO);
	std::replace( date_str.begin(), date_str.end(), ':', '.');
	const std::string name = grid2dpath_in + "/" + date_str + ext;
	openGridFile(name);

	//Radiation parameters
	if(parameter==MeteoGrids::ISWR) readGridLayer("radsw", 2, grid_out);
	if(parameter==MeteoGrids::RSWR) {
		Grid2DObject net;
		readGridLayer("radsw", 2, grid_out);
		readGridLayer("radswnet", 2, net);
		grid_out.grid2D -= net.grid2D;
	}
	if(parameter==MeteoGrids::ILWR) readGridLayer("radlwin", 2, grid_out);
	if(parameter==MeteoGrids::ALB) {
		Grid2DObject rswr, iswr;
		readGridLayer("radsw", 2, iswr);
		readGridLayer("radswnet", 2, grid_out); //net radiation
		rswr.grid2D = iswr.grid2D - grid_out.grid2D;
		grid_out.grid2D = iswr.grid2D/rswr.grid2D;
	}

	//Wind grids
	if(parameter==MeteoGrids::U) readGridLayer("u", 2, grid_out);
	if(parameter==MeteoGrids::V) readGridLayer("v", 2, grid_out);
	if(parameter==MeteoGrids::W) readGridLayer("w", 2, grid_out);
	if(parameter==MeteoGrids::VW) {
		Grid2DObject V;
		readGridLayer("u", 2, grid_out); //U
		readGridLayer("v", 2, V);
		for(unsigned int jj=0; jj<grid_out.nrows; jj++) {
			for(unsigned int ii=0; ii<grid_out.ncols; ii++) {
164
				grid_out(ii,jj) = sqrt( Optim::pow2(grid_out(ii,jj)) + Optim::pow2(V(ii,jj)) );
165
166
167
168
169
170
171
172
173
			}
		}
	}
	if(parameter==MeteoGrids::DW) {
		Grid2DObject V;
		readGridLayer("u", 2, grid_out); //U
		readGridLayer("v", 2, V);
		for(unsigned int jj=0; jj<grid_out.nrows; jj++) {
			for(unsigned int ii=0; ii<grid_out.ncols; ii++) {
174
				grid_out(ii,jj) = fmod( atan2( grid_out(ii,jj), V(ii,jj) ) * Cst::to_deg + 360., 360.); // turn into degrees [0;360)
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
			}
		}
	}

	//Basic meteo parameters
	if(parameter==MeteoGrids::P) readGridLayer("p", 2, grid_out);
	if(parameter==MeteoGrids::TSG) readGridLayer("tsoil", 2, grid_out); //or is it tss for us?
	/*if(parameter==MeteoGrids::RH) {
		//const double epsilon = Cst::gaz_constant_dry_air / Cst::gaz_constant_water_vapor;
		readGridLayer("qv", 2, grid_out); //water vapor mixing ratio
		//Atmosphere::waterSaturationPressure(T);
		//HACK: compute relative humidity out of it!
		//through potential temperature -> local temperature?
	}*/

	//Hydrological parameters
	if(parameter==MeteoGrids::HS) readGridLayer("snowdpth", 2, grid_out);
	if(parameter==MeteoGrids::HNW) {
		readGridLayer("prcrate1", 2, grid_out); //in kg/m^2/s
		grid_out.grid2D *= 3600.; //we need kg/m^2/h
	}

	//DEM
	std::string dem_marker="zp coordinat";
	if(!is_true_arps) dem_marker="zp_coordinat";
	if(parameter==MeteoGrids::DEM) readGridLayer(dem_marker, 2, grid_out);
	if(parameter==MeteoGrids::SLOPE) {
		DEMObject dem;
		dem.setUpdatePpt(DEMObject::SLOPE);
		readGridLayer(dem_marker, 2, dem);
		dem.update();
		grid_out.set(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner, dem.slope);
	}
	if(parameter==MeteoGrids::AZI) {
		DEMObject dem;
		dem.setUpdatePpt(DEMObject::SLOPE);
		readGridLayer(dem_marker, 2, dem);
		dem.update();
		grid_out.set(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner, dem.azi);
	}

	if(grid_out.isEmpty()) {
		stringstream ss;
		ss << "No suitable data found for parameter " << MeteoGrids::getParameterName(parameter) << " ";
		ss << "at time step " << date.toString(Date::ISO) << " in file \"" << name << "\"";
		throw NoAvailableDataException(ss.str(), AT);
	}

	rewind(fin);
224
225
}

226
void ARPSIO::read3DGrid(Grid3DObject& grid_out, const std::string& i_name)
227
{
228
	const std::string _filename = grid2dpath_in + "/" + i_name;
229
230
231
232
233
234
	openGridFile(_filename);

	//resize the grid just in case
	grid_out.grid3D.resize(dimx, dimy, dimz);

	// Read until the parameter is found
235
	std::string parameter; //HACK
236
237
238
239
240
241
242
	moveToMarker(parameter);

	//read the data we are interested in
	for (unsigned int ix = 0; ix < dimx; ix++) {
		for (unsigned int iy = 0; iy < dimy; iy++) {
			for (unsigned int iz = 0; iz < dimz; iz++) {
				double tmp;
243
244
245
246
247
248
				if(fscanf(fin," %16lf%*[\n]",&tmp)==1) {
					grid_out.grid3D(ix,iy,iz) = tmp;
				} else {
					cleanup();
					throw InvalidFormatException("Failure in reading 3D grid in file "+_filename, AT);
				}
249
250
251
252
253
254
255
256
257
258
259
260
261
262
			}
		}
	}

	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

void ARPSIO::readDEM(DEMObject& dem_out)
{
	std::string _filename;
	cfg.getValue("DEMFILE", "Input", _filename);
	openGridFile(_filename);
	if(is_true_arps) {
263
		readGridLayer(std::string("zp coordinat"), 2 ,dem_out);
264
	} else {
265
		readGridLayer(std::string("zp_coordinat"), 2 ,dem_out);
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
	}
}

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

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

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

void ARPSIO::readMeteoData(const Date& /*dateStart*/, const Date& /*dateEnd*/,
Mathias Bavay's avatar
Mathias Bavay committed
288
                           std::vector< std::vector<MeteoData> >& /*vecMeteo*/,
289
                           const size_t&)
290
291
292
293
294
295
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

void ARPSIO::writeMeteoData(const std::vector< std::vector<MeteoData> >& /*vecMeteo*/,
Mathias Bavay's avatar
Mathias Bavay committed
296
                            const std::string&)
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

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

void ARPSIO::write2DGrid(const Grid2DObject& /*grid_in*/, const std::string& /*name*/)
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

314
315
316
317
318
319
void ARPSIO::write2DGrid(const Grid2DObject&, const MeteoGrids::Parameters&, const Date&)
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

320
321
322
323
324
325
void ARPSIO::initializeGRIDARPS()
{
	double v1, v2;

	//go to read the sizes
	moveToMarker("nnx");
326
327
328
329
330
	//finish reading the line and move to the next one
	if(fscanf(fin,"%*[^\n]")!=0) {
		cleanup();
		throw InvalidFormatException("Error in file format of file "+filename, AT);
	}
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
	if (fscanf(fin," %u %u %u \n",&dimx,&dimy,&dimz)!=3) {
		cleanup();
		throw InvalidFormatException("Can not read dimx, dimy, dimz from file "+filename, AT);
	}
	if (dimx==0 || dimy==0 || dimz==0) {
		cleanup();
		throw IndexOutOfBoundsException("Invalid dimx, dimy, dimz from file "+filename, AT);
	}

	//initializing cell size
	moveToMarker("x_coordinate");
	if (fscanf(fin,"%lg %lg",&v1,&v2)!=2) {
		cleanup();
		throw InvalidFormatException("Can not read first two x coordinates from file "+filename, AT);
	}
	const double cellsize_x = v2 - v1;
	moveToMarker("y_coordinate");
	if (fscanf(fin,"%lg %lg",&v1,&v2)!=2) {
		cleanup();
		throw InvalidFormatException("Can not read first two y coordinates from file "+filename, AT);
	}
	const double cellsize_y = v2 - v1;
	if(cellsize_x!=cellsize_y) {
		cleanup();
		throw InvalidFormatException("Only square cells currently supported! Non compliance in file "+filename, AT);
	}
	cellsize = cellsize_y;
358
359
360

	//HACK: zcoords must be read from zp. But they are NOT constant...

361
362
363
364
365
366
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
}

void ARPSIO::initializeTrueARPS(const char curr_line[ARPS_MAX_LINE_LENGTH])
{
	double v1, v2;

	//go to read the sizes
	if (sscanf(curr_line," nx = %u, ny = %u, nz = %u ",&dimx,&dimy,&dimz)!=3) {
		cleanup();
		throw InvalidFormatException("Can not read dimx, dimy, dimz from file "+filename, AT);
	}
	if (dimx==0 || dimy==0 || dimz==0) {
		cleanup();
		throw IndexOutOfBoundsException("Invalid dimx, dimy, dimz from file "+filename, AT);
	}

	//initializing cell size
	moveToMarker("x coordinate");
	if (fscanf(fin,"%lg %lg",&v1,&v2)!=2) {
		cleanup();
		throw InvalidFormatException("Can not read first two x coordinates from file "+filename, AT);
	}
	const double cellsize_x = v2 - v1;
	moveToMarker("y coordinate");
	if (fscanf(fin,"%lg %lg",&v1,&v2)!=2) {
		cleanup();
		throw InvalidFormatException("Can not read first two y coordinates from file "+filename, AT);
	}
	const double cellsize_y = v2 - v1;
	if(cellsize_x!=cellsize_y) {
		cleanup();
		throw InvalidFormatException("Only square cells currently supported! Non compliance in file "+filename, AT);
	}
	cellsize = cellsize_y;
395
396
397
398
399
400
401
402
403
404
405

	moveToMarker("z coordinate");
	while (fscanf(fin,"%lg",&v1)==1) {
		zcoord.push_back( v1 );
	}
	if(zcoord.size()!=dimz) {
		stringstream ss;
		ss << "Expected " << dimz << " z coordinates in file \""+filename+"\", found " << zcoord.size();
		cleanup();
		throw InvalidFormatException(ss.str(), AT);
	}
406
407
}

408
void ARPSIO::openGridFile(const std::string& in_filename)
409
410
{
	unsigned int v1;
411
	filename = in_filename;
412
413
414
415
416
417
418
419
420
421

	if((fin=fopen(filename.c_str(),"r")) == NULL) {
		cleanup();
		throw FileAccessException("Can not open file "+filename, AT);
	}

	//identify if the file is an original arps file or a file modified by ARPSGRID
	char dummy[ARPS_MAX_LINE_LENGTH];
	for (int j=0; j<5; j++) {
		//the first easy difference in the structure happens at line 5
422
423
424
425
		if(fgets(dummy,ARPS_MAX_STRING_LENGTH,fin)==NULL) {
			cleanup();
			throw InvalidFormatException("Fail to read header lines of file "+filename, AT);
		}
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
	}
	if (sscanf(dummy," nx = %u, ny = ", &v1)<1) {
		//this is an ASCII file modified by ARPSGRID
		is_true_arps=false;
		initializeGRIDARPS();
	} else {
		//this is a true ARPS file
		initializeTrueARPS(dummy);
	}

	//come back to the begining of the file
	rewind(fin);
}

void ARPSIO::cleanup() throw()
{
	if (fin!=NULL) {//close fin if open
		fclose(fin);
		fin=NULL;
	}
	/*if (fin.is_open()) {//close fin if open
		fin.close();
	}*/
449
450

	zcoord.clear();
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
}

/** @brief Read a specific layer for a given parameter from the ARPS file
 * @param parameter The parameter to extract. This could be any of the following:
 *        - x_coordinate for getting the X coordinates of the mesh
 *        - y_coordinate for getting the Y coordinates of the mesh
 *        - zp_coordinat for getting the Z coordinates of the mesh
 *        - u for getting the u component of the wind field
 *        - v for getting the v component of the wind field
 *        - w for getting the w component of the wind field
 * @param layer     Index of the layer to extract (1 to dimz)
 * @param grid      [out] grid containing the values. The grid will be resized if necessary.
*/
void ARPSIO::readGridLayer(const std::string& parameter, const unsigned int& layer, Grid2DObject& grid)
{
	if(layer<1 || layer>dimz) {
		cleanup();
		stringstream tmp;
469
		tmp << "Layer " << layer << " does not exist in ARPS file " << filename << " (nr layers=" << dimz << ")";
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
		throw IndexOutOfBoundsException(tmp.str(), AT);
	}

	//resize the grid just in case
	Coords llcorner(coordin, coordinparam);
	llcorner.setXY(xcoord, ycoord, IOUtils::nodata);
	grid.set(dimx, dimy, cellsize, llcorner);

	// Read until the parameter is found
	moveToMarker(parameter);

	// move to the begining of the layer of interest
	if(layer>1) {
		double tmp;
		const unsigned int jmax=dimx*dimy*(layer-1);
		for (unsigned int j = 0; j < jmax; j++)
486
487
488
489
			if(fscanf(fin," %16lf%*[\n]",&tmp)==EOF) {
				cleanup();
				throw InvalidFormatException("Fail to skip data layers in file "+filename, AT);
			}
490
491
492
	}

	//read the data we are interested in
493
494
	for (unsigned int iy = 0; iy < dimy; iy++) {
		for (unsigned int ix = 0; ix < dimx; ix++) {
495
			double tmp;
496
497
498
499
500
501
			if(fscanf(fin," %16lf%*[\n]",&tmp)==1) {
				grid.grid2D(ix,iy) = tmp;
			} else {
				cleanup();
				throw InvalidFormatException("Fail to read data layer in file "+filename, AT);
			}
502
503
504
505
506
507
508
		}
	}
}

void ARPSIO::moveToMarker(const std::string& marker)
{
	char dummy[ARPS_MAX_LINE_LENGTH];
509
	int nb_elems=0;
510
	do {
511
512
		nb_elems=fscanf(fin," %[^\t\n] ",dummy);
	} while (!feof(fin) && strcmp(dummy,marker.c_str()) != 0 && nb_elems!=0);
513
514
515
516
517
	if(feof(fin)) {
		cleanup();
		const std::string message = "End of file "+filename+" should NOT have been reached when looking for "+marker;
		throw InvalidFormatException(message, AT);
	}
518
519
520
521
522
	if(nb_elems==0) {
		cleanup();
		const std::string message = "Matching failure in file "+filename+" when looking for "+marker;
		throw InvalidFormatException(message, AT);
	}
523
524
525
}

} //namespace