WSL/SLF GitLab Repository

Grid2DObject.cc 14.3 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/dataClasses/Grid2DObject.h>
19
20
#include <meteoio/IOExceptions.h>
#include <meteoio/IOUtils.h>
21
22
#include <meteoio/MathOptim.h>
#include <meteoio/ResamplingAlgorithms2D.h>
23
#include <cmath>
24
25
26
27
28
29

using namespace std;

namespace mio {

Grid2DObject& Grid2DObject::operator=(const Grid2DObject& source) {
30
	if (this != &source) {
31
32
33
34
35
36
37
		grid2D = source.grid2D;
		cellsize = source.cellsize;
		llcorner = source.llcorner;
	}
	return *this;
}

38
39
40
41
42
43
44
45
46
47
Grid2DObject& Grid2DObject::operator=(const double& value) {
	grid2D = value;
	return *this;
}

Grid2DObject& Grid2DObject::operator+=(const double& rhs) {
	grid2D += rhs;
	return *this;
}

48
const Grid2DObject Grid2DObject::operator+(const double& rhs) const {
49
50
51
52
53
54
55
56
57
58
59
60
	Grid2DObject result = *this;
	result.grid2D += rhs;
	return result;
}

Grid2DObject& Grid2DObject::operator+=(const Grid2DObject& rhs) {
	if (!isSameGeolocalization(rhs))
		throw InvalidArgumentException("[E] grids must have the same geolocalization in order to do arithmetic operations!", AT);
	grid2D += rhs.grid2D;
	return *this;
}

61
const Grid2DObject Grid2DObject::operator+(const Grid2DObject& rhs) const {
62
63
64
65
66
67
68
69
70
71
72
73
	if (!isSameGeolocalization(rhs))
		throw InvalidArgumentException("[E] grids must have the same geolocalization in order to do arithmetic operations!", AT);
	Grid2DObject result(*this);
	result.grid2D += rhs.grid2D;
	return result;
}

Grid2DObject& Grid2DObject::operator-=(const double& rhs) {
	grid2D -= rhs;
	return *this;
}

74
const Grid2DObject Grid2DObject::operator-(const double& rhs) const {
75
76
77
78
79
80
81
82
83
84
85
86
	Grid2DObject result(*this);
	result.grid2D -= rhs;
	return result;
}

Grid2DObject& Grid2DObject::operator-=(const Grid2DObject& rhs) {
	if (!isSameGeolocalization(rhs))
		throw InvalidArgumentException("[E] grids must have the same geolocalization in order to do arithmetic operations!", AT);
	grid2D -= rhs.grid2D;
	return *this;
}

87
const Grid2DObject Grid2DObject::operator-(const Grid2DObject& rhs) const {
88
89
90
91
92
93
94
95
96
97
98
99
	if (!isSameGeolocalization(rhs))
		throw InvalidArgumentException("[E] grids must have the same geolocalization in order to do arithmetic operations!", AT);
	Grid2DObject result(*this);
	result.grid2D -= rhs.grid2D;
	return result;
}

Grid2DObject& Grid2DObject::operator*=(const double& rhs) {
	grid2D *= rhs;
	return *this;
}

100
const Grid2DObject Grid2DObject::operator*(const double& rhs) const {
101
102
103
104
105
106
107
108
109
110
111
112
	Grid2DObject result(*this);
	result.grid2D *= rhs;
	return result;
}

Grid2DObject& Grid2DObject::operator*=(const Grid2DObject& rhs) {
	if (!isSameGeolocalization(rhs))
		throw InvalidArgumentException("[E] grids must have the same geolocalization in order to do arithmetic operations!", AT);
	grid2D *= rhs.grid2D;
	return *this;
}

113
const Grid2DObject Grid2DObject::operator*(const Grid2DObject& rhs) const {
114
115
116
117
118
119
120
121
122
123
124
125
	if (!isSameGeolocalization(rhs))
		throw InvalidArgumentException("[E] grids must have the same geolocalization in order to do arithmetic operations!", AT);
	Grid2DObject result(*this);
	result.grid2D *= rhs.grid2D;
	return result;
}

Grid2DObject& Grid2DObject::operator/=(const double& rhs) {
	grid2D /= rhs;
	return *this;
}

126
const Grid2DObject Grid2DObject::operator/(const double& rhs) const {
127
128
129
130
131
132
133
134
135
136
137
138
	Grid2DObject result(*this);
	result.grid2D /= rhs;
	return result;
}

Grid2DObject& Grid2DObject::operator/=(const Grid2DObject& rhs) {
	if (!isSameGeolocalization(rhs))
		throw InvalidArgumentException("[E] grids must have the same geolocalization in order to do arithmetic operations!", AT);
	grid2D /= rhs.grid2D;
	return *this;
}

139
const Grid2DObject Grid2DObject::operator/(const Grid2DObject& rhs) const {
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
	if (!isSameGeolocalization(rhs))
		throw InvalidArgumentException("[E] grids must have the same geolocalization in order to do arithmetic operations!", AT);
	Grid2DObject result(*this);
	result.grid2D /= rhs.grid2D;
	return result;
}

bool Grid2DObject::operator==(const Grid2DObject& in) const {
	return (isSameGeolocalization(in) && grid2D==in.grid2D);
}

bool Grid2DObject::operator!=(const Grid2DObject& in) const {
	return !(*this==in);
}

155
156
157
158
/*
 * Default constructor.
 * grid2D attribute is initialized by Array2D default constructor.
 */
159
Grid2DObject::Grid2DObject() : grid2D(), llcorner(), cellsize(0.) {}
160

161
Grid2DObject::Grid2DObject(const size_t& i_ncols, const size_t& i_nrows,
162
163
                           const double& i_cellsize, const Coords& i_llcorner)
              : grid2D(i_ncols, i_nrows, IOUtils::nodata), llcorner(i_llcorner),
164
                cellsize(i_cellsize) {}
165

166
167
Grid2DObject::Grid2DObject(const double& i_cellsize, const Coords& i_llcorner, const Array2D<double>& i_grid2D)
              : grid2D(i_grid2D), llcorner(i_llcorner), cellsize(i_cellsize) {}
168

169
Grid2DObject::Grid2DObject(const size_t& i_ncols, const size_t& i_nrows,
170
                           const double& i_cellsize, const Coords& i_llcorner, const double& init)
171
              : grid2D(i_ncols, i_nrows, init), llcorner(i_llcorner), cellsize(i_cellsize) {}
172
173

Grid2DObject::Grid2DObject(const Grid2DObject& i_grid, const double& init)
174
175
              : grid2D(i_grid.grid2D.getNx(), i_grid.grid2D.getNy(), init),
                llcorner(i_grid.llcorner), cellsize(i_grid.cellsize) {}
176

177
178
Grid2DObject::Grid2DObject(const Grid2DObject& i_grid2Dobj, const size_t& i_nx, const size_t& i_ny,
                           const size_t& i_ncols, const size_t& i_nrows)
179
              : grid2D(i_grid2Dobj.grid2D, i_nx,i_ny, i_ncols,i_nrows), llcorner(i_grid2Dobj.llcorner), cellsize(i_grid2Dobj.cellsize)
180
181
182
{
	//we take the previous corner (so we use the same projection parameters)
	//and we shift it by the correct X and Y distance
183
	//llcorner = i_grid2Dobj.llcorner;
184
	if ( (llcorner.getEasting()!=IOUtils::nodata) && (llcorner.getNorthing()!=IOUtils::nodata) ) {
Mathias Bavay's avatar
Mathias Bavay committed
185
186
		llcorner.setXY( llcorner.getEasting()+static_cast<double>(i_nx)*i_grid2Dobj.cellsize,
		                llcorner.getNorthing()+static_cast<double>(i_ny)*i_grid2Dobj.cellsize, IOUtils::nodata);
187
188
189
	}
}

190
bool Grid2DObject::gridify(std::vector<Coords>& vec_points, const bool& keep_invalid) const 
191
{
192
193
194
	bool status=true;
	std::vector<Coords>::iterator v_Itr = vec_points.begin();
	while ( v_Itr != vec_points.end() ) {
195
		if ( gridify(*v_Itr)==false ) {
196
			if (!keep_invalid) v_Itr = vec_points.erase(v_Itr);
197
			else ++v_Itr;
198
199
			status=false;
		} else {
200
			++v_Itr;
201
202
203
204
205
206
207
208
209
		}
	}

	return status;
}

bool Grid2DObject::gridify(Coords& point) const {
	std::string proj_type, proj_args;
	point.getProj(proj_type, proj_args);
210
	if (proj_type=="NULL") {
211
212
213
214
		//if the projection was "NULL", we set it to the grid's
		point.copyProj(llcorner);
	}

215
	if (point.getGridI()!=IOUtils::inodata && point.getGridJ()!=IOUtils::inodata) {
216
		//we need to compute (easting,northing) and (lat,lon)
217
218
219
220
221
222
223
224
		return( grid_to_WGS84(point) );
	} else {
		//we need to compute (i,j)
		return( WGS84_to_grid(point) );
	}
}

bool Grid2DObject::grid_to_WGS84(Coords& point) const {
225
	int i = point.getGridI(), j = point.getGridJ();
226

227
	if (i==IOUtils::inodata || j==IOUtils::inodata) {
228
		//the point is invalid (outside the grid or contains nodata)
229
		point.setGridIndex(IOUtils::inodata, IOUtils::inodata, IOUtils::inodata, false); //mark the point as invalid
230
231
232
		return false;
	}

233
234
	const int ncols = (signed)getNx();
	const int nrows = (signed)getNy();
235
	if (i>ncols || i<0 || j>nrows || j<0) {
236
237
		//the point is outside the grid, we reset the indices to the closest values
		//still fitting in the grid and return an error
238
239
240
241
242
		if (i<0) i=0;
		if (j<0) j=0;
		if (i>ncols) i=ncols;
		if (j>nrows) j=nrows;
		point.setGridIndex(i, j, IOUtils::inodata, false); //mark the point as invalid
243
244
245
246
247
248
249
		return false;
	}

	//easting and northing in the grid's projection
	const double easting = ((double)i) * cellsize + llcorner.getEasting();
	const double northing = ((double)j) * cellsize + llcorner.getNorthing();

250
	if (point.isSameProj(llcorner)==true) {
251
252
253
254
255
256
257
258
259
260
		//same projection between the grid and the point -> precise, simple and efficient arithmetics
		point.setXY(easting, northing, IOUtils::nodata);
	} else {
		//projections are different, so we have to do an intermediate step...
		Coords tmp_proj;
		tmp_proj.copyProj(point); //making a copy of the original projection
		point.copyProj(llcorner); //taking the grid's projection
		point.setXY(easting, northing, IOUtils::nodata);
		point.copyProj(tmp_proj); //back to the original projection -> reproject the coordinates
	}
261

262
	point.setGridIndex(i, j, IOUtils::inodata, true); //mark the point as valid
263
264
265
266
	return true;
}

bool Grid2DObject::WGS84_to_grid(Coords& point) const {
267
	if (point.getLat()==IOUtils::nodata || point.getLon()==IOUtils::nodata) {
268
			//if the point is invalid, there is nothing we can do
269
			point.setGridIndex(IOUtils::inodata, IOUtils::inodata, IOUtils::inodata, false); //mark the point as invalid
270
271
272
			return false;
	}

273
	bool status=true;
274
275
	int i,j;

276
	if (point.isSameProj(llcorner)==true) {
277
278
279
280
281
282
283
284
285
286
287
288
289
		//same projection between the grid and the point -> precise, simple and efficient arithmetics
		i = (int)floor( (point.getEasting()-llcorner.getEasting()) / cellsize );
		j = (int)floor( (point.getNorthing()-llcorner.getNorthing()) / cellsize );
	} else {
		//projections are different, so we have to do an intermediate step...
		Coords tmp_point(point);
		tmp_point.copyProj(llcorner); //getting the east/north coordinates in the grid's projection
		i = (int)floor( (tmp_point.getEasting()-llcorner.getEasting()) / cellsize );
		j = (int)floor( (tmp_point.getNorthing()-llcorner.getNorthing()) / cellsize );
	}

	//checking that the calculated indices fit in the grid2D
	//and giving them the closest value within the grid if not.
290
	if (i<0) {
291
		i=0;
292
		status=false;
293
	}
294
	if (i>=(signed)getNx()) {
295
		i=(signed)getNx();
296
		status=false;
297
	}
298
	if (j<0) {
299
		j=0;
300
		status=false;
301
	}
302
	if (j>=(signed)getNy()) {
303
		j=(signed)getNy();
304
		status=false;
305
306
	}

307
308
	point.setGridIndex(i, j, IOUtils::inodata, status); //mark as valid or invalid according to status
	return status;
309
310
}

311
void Grid2DObject::set(const size_t& i_ncols, const size_t& i_nrows,
312
                       const double& i_cellsize, const Coords& i_llcorner)
313
{
314
	grid2D.resize(i_ncols, i_nrows, IOUtils::nodata);
315
	setValues(i_cellsize, i_llcorner);
316
317
}

318
void Grid2DObject::set(const size_t& i_ncols, const size_t& i_nrows,
319
320
                       const double& i_cellsize, const Coords& i_llcorner, const double& init)
{
321
	grid2D.resize(i_ncols, i_nrows, init);
322
	setValues(i_cellsize, i_llcorner);
323
324
}

325
void Grid2DObject::set(const double& i_cellsize, const Coords& i_llcorner, const Array2D<double>& i_grid2D)
326
{
327
	setValues(i_cellsize, i_llcorner);
328
	grid2D = i_grid2D;
329
330
}

331
332
333
334
335
336
void Grid2DObject::set(const Grid2DObject& i_grid, const double& init)
{
	setValues(i_grid.cellsize, i_grid.llcorner);
	grid2D.resize(i_grid.grid2D.getNx(), i_grid.grid2D.getNy(), init);
}

337
338
339
340
341
342
343
344
345
346
347
void Grid2DObject::rescale(const double& i_cellsize)
{
	if (grid2D.getNx()==0 || grid2D.getNy()==0)
		throw InvalidArgumentException("Can not rescale an empty grid!", AT);
	
	const double factor_x = cellsize / i_cellsize;
	const double factor_y = cellsize / i_cellsize;
	grid2D = ResamplingAlgorithms2D::BilinearResampling(grid2D, factor_x, factor_y);
	cellsize = i_cellsize;
}

348
void Grid2DObject::size(size_t& o_ncols, size_t& o_nrows) const {
349
350
	o_ncols = getNx();
	o_nrows = getNy();
351
352
}

353
354
355
356
size_t Grid2DObject::size() const {
	return grid2D.size();
}

357
size_t Grid2DObject::getNx() const {
358
	return grid2D.getNx();
359
360
}

361
size_t Grid2DObject::getNy() const {
362
	return grid2D.getNy();
363
364
365
}


366
367
368
369
void Grid2DObject::clear() {
	grid2D.clear();
}

370
bool Grid2DObject::empty() const {
371
	return (grid2D.getNx()==0 && grid2D.getNy()==0);
372
373
}

374
void Grid2DObject::setValues(const double& i_cellsize, const Coords& i_llcorner)
375
{
376
377
	cellsize = i_cellsize;
	llcorner = i_llcorner;
378
379
}

380
bool Grid2DObject::isSameGeolocalization(const Grid2DObject& target) const
381
{
382
	if ( grid2D.getNx()==target.grid2D.getNx() &&
383
384
385
	    grid2D.getNy()==target.grid2D.getNy() &&
	    llcorner==target.llcorner &&
	    cellsize==target.cellsize) {
386
387
388
389
390
391
		return true;
	} else {
		return false;
	}
}

392
393
bool Grid2DObject::clusterization(const std::vector<double>& thresholds, const std::vector<double>& ids)
{
394
	if (thresholds.empty()==0) {
395
396
397
398
399
		throw IOException("Can't start clusterization, cluster definition list is empty", AT);
	}
	if ((thresholds.size()+1) != ids.size()) {
		throw IOException("Can't start clusterization, cluster definition list doesnt fit id definition list", AT);
	}
400
	const size_t nscl = thresholds.size();
401
	for (size_t jj = 0; jj< grid2D.size(); jj++){
402
403
		const double& val = grid2D(jj);
		if (val!=IOUtils::nodata){
404
			size_t i = 0;
405
			for ( ;i<nscl; i++)
406
				if (thresholds[i] >= val)
407
408
409
410
411
412
413
					break;
			grid2D(jj) = ids[i];
		}
	}
	return true;
}

414
double& Grid2DObject::operator()(const size_t& ix, const size_t& iy) {
415
416
	return grid2D(ix,iy);
}
417

418
double Grid2DObject::operator()(const size_t& ix, const size_t& iy) const {
419
420
	return grid2D(ix,iy);
}
421

422
double& Grid2DObject::operator()(const size_t& i) {
423
424
425
	return grid2D(i);
}

426
double Grid2DObject::operator()(const size_t& i) const {
427
428
429
	return grid2D(i);
}

430
const std::string Grid2DObject::toString() const {
431
	std::ostringstream os;
432
	os << "<Grid2DObject>\n";
433
	os << llcorner.toString();
434
	os << grid2D.getNx() << " x " << grid2D.getNy() << " @ " << cellsize << "m\n";
435
436
437
438
439
440
441
	os << grid2D.toString();
	os << "</Grid2DObject>\n";
	return os.str();
}

std::iostream& operator<<(std::iostream& os, const Grid2DObject& grid) {
	os.write(reinterpret_cast<const char*>(&grid.cellsize), sizeof(grid.cellsize));
442
443
444
445
446
	os << grid.llcorner;
	os << grid.grid2D;
	return os;
}

447
448
449
450
std::iostream& operator>>(std::iostream& is, Grid2DObject& grid) {
	is.read(reinterpret_cast<char*>(&grid.cellsize), sizeof(grid.cellsize));
	is >> grid.llcorner;
	is >> grid.grid2D;
451
	return is;
452
453
}

454
} //namespace