WSL/SLF GitLab Repository

DEMObject.h 6.76 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
/***********************************************************************************/
/*  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/>.
*/
#ifndef DEMOBJECT_H
#define DEMOBJECT_H

#include <meteoio/Array2D.h>
#include <meteoio/Grid2DObject.h>
#include <meteoio/IOUtils.h>
24
#include <meteoio/meteolaws/Meteoconst.h> //for math constants
25

26
#include <cmath>
27
28
29
30
31
32
#include <limits>

namespace mio {

/**
 * @class DEMObject
33
34
35
36
 * @brief A class to represent DEMs and automatically compute some properties.
 * This class stores elevation grids and their georeferencing, expressed as the lower-left coordinates, the cellsize (cells are assumed to be square) and a nodata code for potentially empty cells (The nodata parameter is supposed to be IOUtils::nodata).
 * This class also automatically computes local slope, azimuth, curvature, normals and minimal/maximal for normalization.
 * Various algorithms are available to compute these properties (see mio::DEMObject::slope_type) and it is possible to toggle between automatic refresh or not. Several other DEM related values can be computed, such as the horizon, displacements within the DEM, etc
37
38
 *
 * @ingroup data_str
39
40
41
 * @author Gaël Rosset - Mathias Bavay
 * @date   2009-07-20
 */
42
43
44
45
46
47
48
49
50

#ifdef _POPC_
class DEMObjectDummy {}; //HACK for POPC

#include <paroc_base.h>
class DEMObject : public Grid2DObject/*, POPBase*/ {
	public:
		void Serialize(POPBuffer &buf, bool pack);
#else
51
class DEMObject : public Grid2DObject {
52
#endif
53
54
55
56
57
	public:
		Array2D<double> slope;
		Array2D<double> azi;
		Array2D<double> curvature;
		Array2D<double> Nx, Ny, Nz;
58
59
		double min_altitude, min_slope, min_curvature;
		double max_altitude, max_slope, max_curvature;
60
61
62
63
64
65
66

		///Keywords for slope computation algorithm
		typedef enum SLOPE_TYPE {
			DFLT, ///< whatever algorithm that has been defined as default
			FLEM, ///< four nearest neighbors (Fleming and Hoffer, 1979). It seems to be the same as (Zevenbergen and Thorne, 1987)
			HICK, ///< maximum downhill slope method (Dunn and Hickey, 1998)
			HORN, ///< eight neighbor algorithm (Horn, 1981) as used by ArcGIS. It seems to be the same as (Corripio, 2002) but for border cells.
67
			CORR, ///< surface normal vector using the two triangle method (Corripio, 2003) and eight-neighbor algorithm (Horn, 1981) for border cells
68
69
			D8 ///< discretized azimuth directions (angles for N, NE, etc) and slope rounded to nearest integer
		} slope_type;
70
71
72
73
74
75
76
77

		///Keywords for automatic update of parameters. They can be combined with "|"
		typedef enum UPDATE_TYPE {
			NO_UPDATE=0, ///< no updates at all
			SLOPE=1, ///< update the slopes
			NORMAL=2, ///< update the normals
			CURVATURE=4 ///< update the curvatures
		} update_type;
78

79
		DEMObject(const slope_type& i_algorithm=DFLT);
80

Mathias Bavay's avatar
Mathias Bavay committed
81
		DEMObject(const size_t& ncols_in, const size_t& nrows_in,
82
		          const double& cellsize_in, const Coords& llcorner_in, const slope_type& i_algorithm=DFLT);
83

Mathias Bavay's avatar
Mathias Bavay committed
84
		DEMObject(const size_t& ncols_in, const size_t& nrows_in,
85
		          const double& cellsize_in, const Coords& llcorner_in, const Array2D<double>& altitude_in,
86
		          const bool& i_update=true, const slope_type& i_algorithm=DFLT);
87

88
		DEMObject(const Grid2DObject& dem_in, const bool& i_update=true, const slope_type& i_algorithm=DFLT);
89

90
		DEMObject (const DEMObject& i_dem,
Mathias Bavay's avatar
Mathias Bavay committed
91
92
		           const size_t& i_nx, const size_t& i_ny, //Point in the plane
		           const size_t& i_ncols, const size_t& i_nrows, //dimensions of the sub-plane
93
		           const bool& i_update=true, const slope_type& i_algorithm=DFLT);
94

95
		void setDefaultAlgorithm(const slope_type& i_algorithm);
96
		int getDefaultAlgorithm() const;
97
		void setUpdatePpt(const update_type& in_update_flag);
98
		int getUpdatePpt() const;
99

100
101
102
103
104
105
106
107
108
		void update(const std::string& algorithm);
		void update(const slope_type& algorithm=DFLT);
		void updateAllMinMax();
		void printFailures();
		void sanitize();
		double horizontalDistance(const double& xcoord1, const double& ycoord1, const double& xcoord2, const double& ycoord2);
		double horizontalDistance(Coords point1, const Coords& point2);
		double terrainDistance(Coords point1, const Coords& point2);
		void getPointsBetween(Coords point1, Coords point2, std::vector<GRID_POINT_2D>& vec_points);
109
		void getPointsBetween(const Coords& point, const double& bearing, std::vector<GRID_POINT_2D>& vec_points);
110
111
		double getHorizon(const Coords& point, const double& bearing);
		void getHorizon(const Coords& point, const double& increment, std::vector<double>& horizon);
112

113
114
115
		friend std::iostream& operator<<(std::iostream& os, const DEMObject& dem);
		friend std::iostream& operator>>(std::iostream& is, DEMObject& dem);

116
117
	private:
		void CalculateAziSlopeCurve(slope_type algorithm);
118
119
120
121
122
123
		double CalculateAspect(const double& o_Nx, const double& o_Ny, const double& o_Nz, const double& o_slope, const double no_slope=Cst::PI);
		void CalculateHick(double A[4][4], double& o_slope, double& o_Nx, double& o_Ny, double& o_Nz);
		void CalculateFleming(double A[4][4], double& o_slope, double& o_Nx, double& o_Ny, double& o_Nz);
		void CalculateHorn(double A[4][4], double& o_slope, double& o_Nx, double& o_Ny, double& o_Nz);
		void CalculateCorripio(double A[4][4], double& o_slope, double& o_Nx, double& o_Ny, double& o_Nz);
		void (DEMObject::*CalculateSlope)(double A[4][4], double& o_slope, double& o_Nx, double& o_Ny, double& o_Nz);
124
125
126
127
128
129
130
		double getCurvature(double A[4][4]);

		double steepestGradient(double A[4][4]);
		double lineGradient(const double& A1, const double& A2, const double& A3);
		double fillMissingGradient(const double& delta1, const double& delta2);
		void surfaceGradient(double& dx_sum, double& dy_sum, double A[4][4]);
		double avgHeight(const double& z1, const double &z2, const double& z3);
Mathias Bavay's avatar
Mathias Bavay committed
131
		void getNeighbours(const size_t i, const size_t j, double A[4][4]);
132
133
		double safeGet(const int i, const int j);

134
		int update_flag;
135
		slope_type dflt_algorithm;
Mathias Bavay's avatar
Mathias Bavay committed
136
137
		size_t slope_failures; ///<contains the number of points that have an elevation but no slope
		size_t curvature_failures; ///<contains the number of points that have an elevation but no curvature
138
139
140
141
};
} //end namespace

#endif