WSL/SLF GitLab Repository

DEMObject.h 6.29 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
33
34
#include <limits>

namespace mio {

/**
 * @class DEMObject
 * @brief A class to represent DEMs: reads elevation grids, computes local slope, azimuth, curvature.
 * The nodata parameter is supposed to be IOUtils::nodata.
35
36
 *
 * @ingroup data_str
37
38
39
 * @author Gaël Rosset - Mathias Bavay
 * @date   2009-07-20
 */
40
41
42
43
44
45
46
47
48

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

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

		///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.
65
			CORR, ///< surface normal vector using the two triangle method (Corripio, 2003) and eight-neighbor algorithm (Horn, 1981) for border cells
66
67
			D8 ///< discretized azimuth directions (angles for N, NE, etc) and slope rounded to nearest integer
		} slope_type;
68
69
70
71
72
73
74
75

		///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;
76

77
		DEMObject(const slope_type& i_algorithm=DFLT);
78

79
		DEMObject(const unsigned int& ncols_in, const unsigned int& nrows_in,
80
		          const double& cellsize_in, const Coords& llcorner_in, const slope_type& i_algorithm=DFLT);
81

82
		DEMObject(const unsigned int& ncols_in, const unsigned int& nrows_in,
83
		          const double& cellsize_in, const Coords& llcorner_in, const Array2D<double>& altitude_in,
84
		          const bool& i_update=true, const slope_type& i_algorithm=DFLT);
85

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

88
89
90
91
		DEMObject (const DEMObject& i_dem,
		           const unsigned int& i_nx, const unsigned int& i_ny, //Point in the plane
		           const unsigned int& i_ncols, const unsigned int& i_nrows, //dimensions of the sub-plane
		           const bool& i_update=true, const slope_type& i_algorithm=DFLT);
92

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

98
99
100
101
102
103
104
105
106
		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);
107
		void getPointsBetween(const Coords& point, const double& bearing, std::vector<GRID_POINT_2D>& vec_points);
108
109
		double getHorizon(const Coords& point, const double& bearing);
		void getHorizon(const Coords& point, const double& increment, std::vector<double>& horizon);
110

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

114
115
	private:
		void CalculateAziSlopeCurve(slope_type algorithm);
116
117
118
119
120
121
		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);
122
123
124
125
126
127
128
129
130
131
		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);
		void getNeighbours(const unsigned int i, const unsigned int j, double A[4][4]);
		double safeGet(const int i, const int j);

132
		int update_flag;
133
		slope_type dflt_algorithm;
134
135
136
137
138
139
		unsigned int slope_failures; ///<contains the number of points that have an elevation but no slope
		unsigned int curvature_failures; ///<contains the number of points that have an elevation but no curvature
};
} //end namespace

#endif