WSL/SLF GitLab Repository

Commit a97ddf12 authored by Mathias Bavay's avatar Mathias Bavay
Browse files

Better UTM zone check, some fixes for variables shadowing and alignment

parent 13a96861
......@@ -677,11 +677,11 @@ short int Coords::getEPSG() const {
char zoneLetter;
parseUTMZone(coordparam, zoneLetter, zoneNumber);
if(zoneLetter >= 'M') {
//northern hemisphere
return (32600+zoneNumber);
//northern hemisphere. We KNOW it will fit in a short int
return (static_cast<short int>(32600+zoneNumber));
} else {
//southern hemisphere
return (32700+zoneNumber);
//southern hemisphere. We KNOW it will fit in a short int
return (static_cast<short int>(32700+zoneNumber));
}
}
if(coordsystem=="UPS") {
......@@ -1378,6 +1378,9 @@ void Coords::parseUTMZone(const std::string& zone_info, char& zoneLetter, short
(sscanf(zone_info.c_str(), "%hd %c)", &zoneNumber, &zoneLetter) < 2)) {
throw InvalidFormatException("Can not parse given UTM zone: "+zone_info,AT);
}
if(zoneNumber<1 || zoneNumber>60) {
throw InvalidFormatException("Invalid UTM zone: "+zone_info+" (zone should be between 1 and 60, zone letter in [C-N, P-X])",AT);
}
zoneLetter = (char)toupper(zoneLetter); //just in case... (sorry for the pun!)
if(zoneLetter=='Y' || zoneLetter=='Z' || zoneLetter=='A' || zoneLetter=='B') {
//Special zones for the poles: we should NOT use UTM in these regions!
......
......@@ -708,29 +708,29 @@ void DEMObject::CalculateAziSlopeCurve(slope_type algorithm) {
} // end of CalculateAziSlope
double DEMObject::CalculateAspect(const double& Nx, const double& Ny, const double& Nz, const double& slope, const double no_slope) {
double DEMObject::CalculateAspect(const double& o_Nx, const double& o_Ny, const double& o_Nz, const double& o_slope, const double no_slope) {
//Calculates the aspect at a given point knowing its normal vector and slope
//(direction of the normal pointing out of the surface, clockwise from north)
//This azimuth calculation is similar to Hodgson (1998)
//local_nodata is the value that we want to give to the aspect of points that don't have a slope
//The value is a bearing (ie: deg, clockwise, 0=North)
if(Nx==IOUtils::nodata || Ny==IOUtils::nodata || Nz==IOUtils::nodata || slope==IOUtils::nodata) {
if(o_Nx==IOUtils::nodata || o_Ny==IOUtils::nodata || o_Nz==IOUtils::nodata || o_slope==IOUtils::nodata) {
return IOUtils::nodata;
}
if ( slope > 0. ) { //there is some slope
if ( Nx == 0. ) { //no E-W slope, so it is purely N-S
if ( Ny < 0. ) {
if ( o_slope > 0. ) { //there is some slope
if ( o_Nx == 0. ) { //no E-W slope, so it is purely N-S
if ( o_Ny < 0. ) {
return(180.); // south facing
} else {
return (0.); // north facing
}
} else { //there is a E-W slope
if ( Nx > 0. ) {
return (90. - atan(Ny/Nx)*Cst::to_deg);
if ( o_Nx > 0. ) {
return (90. - atan(o_Ny/o_Nx)*Cst::to_deg);
} else {
return (270. - atan(Ny/Nx)*Cst::to_deg);
return (270. - atan(o_Ny/o_Nx)*Cst::to_deg);
}
}
} else { // if slope = 0
......@@ -739,88 +739,88 @@ double DEMObject::CalculateAspect(const double& Nx, const double& Ny, const doub
}
void DEMObject::CalculateHick(double A[4][4], double& slope, double& Nx, double& Ny, double& Nz) {
void DEMObject::CalculateHick(double A[4][4], double& o_slope, double& o_Nx, double& o_Ny, double& o_Nz) {
//This calculates the surface normal vector using the steepest slope method (Dunn and Hickey, 1998):
//the steepest slope found in the eight cells surrounding (i,j) is given to be the slope in (i,j)
//Beware, sudden steps could happen
const double smax = steepestGradient(A); //steepest local gradient
if(smax==IOUtils::nodata) {
slope = IOUtils::nodata;
Nx = IOUtils::nodata;
Ny = IOUtils::nodata;
Nz = IOUtils::nodata;
o_slope = IOUtils::nodata;
o_Nx = IOUtils::nodata;
o_Ny = IOUtils::nodata;
o_Nz = IOUtils::nodata;
slope_failures++;
} else {
slope = atan(smax)*Cst::to_deg;
o_slope = atan(smax)*Cst::to_deg;
//Nx and Ny: x and y components of the normal pointing OUT of the surface
if ( smax > 0. ) { //ie: there is some slope
double dx_sum, dy_sum;
surfaceGradient(dx_sum, dy_sum, A);
if(dx_sum==IOUtils::nodata || dy_sum==IOUtils::nodata) {
Nx = IOUtils::nodata;
Ny = IOUtils::nodata;
Nz = IOUtils::nodata;
o_Nx = IOUtils::nodata;
o_Ny = IOUtils::nodata;
o_Nz = IOUtils::nodata;
slope_failures++;
} else {
Nx = -1.0 * dx_sum / (2. * cellsize); //Nx=-dz/dx
Ny = -1.0 * dy_sum / (2. * cellsize); //Ny=-dz/dy
Nz = 1.; //Nz=1 (normalized by definition of Nx and Ny)
o_Nx = -1.0 * dx_sum / (2. * cellsize); //Nx=-dz/dx
o_Ny = -1.0 * dy_sum / (2. * cellsize); //Ny=-dz/dy
o_Nz = 1.; //Nz=1 (normalized by definition of Nx and Ny)
}
} else { //ie: there is no slope
Nx = 0.;
Ny = 0.;
Nz = 1.;
o_Nx = 0.;
o_Ny = 0.;
o_Nz = 1.;
}
}
}
void DEMObject::CalculateFleming(double A[4][4], double& slope, double& Nx, double& Ny, double& Nz) {
void DEMObject::CalculateFleming(double A[4][4], double& o_slope, double& o_Nx, double& o_Ny, double& o_Nz) {
//This calculates the surface normal vector using method by Fleming and Hoffer (1979)
if(A[2][1]!=IOUtils::nodata && A[2][3]!=IOUtils::nodata && A[3][2]!=IOUtils::nodata && A[1][2]!=IOUtils::nodata) {
Nx = 0.5 * (A[2][1] - A[2][3]) / cellsize;
Ny = 0.5 * (A[3][2] - A[1][2]) / cellsize;
Nz = 1.;
slope = atan( sqrt(Nx*Nx+Ny*Ny) ) * Cst::to_deg;
o_Nx = 0.5 * (A[2][1] - A[2][3]) / cellsize;
o_Ny = 0.5 * (A[3][2] - A[1][2]) / cellsize;
o_Nz = 1.;
o_slope = atan( sqrt(o_Nx*o_Nx+o_Ny*o_Ny) ) * Cst::to_deg;
} else {
CalculateHick(A, slope, Nx, Ny, Nz);
CalculateHick(A, o_slope, o_Nx, o_Ny, o_Nz);
}
}
void DEMObject::CalculateHorn(double A[4][4], double& slope, double& Nx, double& Ny, double& Nz) {
void DEMObject::CalculateHorn(double A[4][4], double& o_slope, double& o_Nx, double& o_Ny, double& o_Nz) {
//This calculates the slope using the two eight neighbors method given in Horn (1981)
//This is also the algorithm used by ArcGIS
if ( A[1][1]!=IOUtils::nodata && A[1][2]!=IOUtils::nodata && A[1][3]!=IOUtils::nodata &&
A[2][1]!=IOUtils::nodata && A[2][2]!=IOUtils::nodata && A[2][3]!=IOUtils::nodata &&
A[3][1]!=IOUtils::nodata && A[3][2]!=IOUtils::nodata && A[3][3]!=IOUtils::nodata) {
Nx = ((A[3][3]+2.*A[2][3]+A[1][3]) - (A[3][1]+2.*A[2][1]+A[1][1])) / (8.*cellsize);
Ny = ((A[1][3]+2.*A[1][2]+A[1][1]) - (A[3][3]+2.*A[3][2]+A[3][1])) / (8.*cellsize);
Nz = 1.;
o_Nx = ((A[3][3]+2.*A[2][3]+A[1][3]) - (A[3][1]+2.*A[2][1]+A[1][1])) / (8.*cellsize);
o_Ny = ((A[1][3]+2.*A[1][2]+A[1][1]) - (A[3][3]+2.*A[3][2]+A[3][1])) / (8.*cellsize);
o_Nz = 1.;
//There is no difference between slope = acos(n_z/|n|) and slope = atan(sqrt(sx*sx+sy*sy))
//slope = acos( (Nz / sqrt( Nx*Nx + Ny*Ny + Nz*Nz )) );
slope = atan( sqrt(Nx*Nx+Ny*Ny) ) * Cst::to_deg;
o_slope = atan( sqrt(o_Nx*o_Nx+o_Ny*o_Ny) ) * Cst::to_deg;
} else {
//steepest slope method (Dunn and Hickey, 1998)
CalculateHick(A, slope, Nx, Ny, Nz);
CalculateHick(A, o_slope, o_Nx, o_Ny, o_Nz);
}
}
void DEMObject::CalculateCorripio(double A[4][4], double& slope, double& Nx, double& Ny, double& Nz) {
void DEMObject::CalculateCorripio(double A[4][4], double& o_slope, double& o_Nx, double& o_Ny, double& o_Nz) {
//This calculates the surface normal vector using the two triangle method given in Corripio (2002)
if ( A[1][2]!=IOUtils::nodata && A[1][3]!=IOUtils::nodata && A[2][2]!=IOUtils::nodata && A[2][3]!=IOUtils::nodata) {
// See Corripio (2002), knowing that here we normalize the result (divided by Nz=cellsize*cellsize)
Nx = (0.5 * (A[2][2] - A[2][3] + A[1][2] - A[1][3]) ) / cellsize;
Ny = (0.5 * (A[2][2] + A[2][3] - A[1][2] - A[1][3]) ) / cellsize;
Nz = 1.;
o_Nx = (0.5 * (A[2][2] - A[2][3] + A[1][2] - A[1][3]) ) / cellsize;
o_Ny = (0.5 * (A[2][2] + A[2][3] - A[1][2] - A[1][3]) ) / cellsize;
o_Nz = 1.;
//There is no difference between slope = acos(n_z/|n|) and slope = atan(sqrt(sx*sx+sy*sy))
//slope = acos( (Nz / sqrt( Nx*Nx + Ny*Ny + Nz*Nz )) );
slope = atan( sqrt(Nx*Nx+Ny*Ny) ) * Cst::to_deg;
o_slope = atan( sqrt(o_Nx*o_Nx+o_Ny*o_Ny) ) * Cst::to_deg;
} else {
//steepest slope method (Dunn and Hickey, 1998)
CalculateHick(A, slope, Nx, Ny, Nz);
CalculateHick(A, o_slope, o_Nx, o_Ny, o_Nz);
}
}
......
......@@ -103,12 +103,12 @@ class DEMObject : public Grid2DObject {
private:
void CalculateAziSlopeCurve(slope_type algorithm);
double CalculateAspect(const double& Nx, const double& Ny, const double& Nz, const double& slope, const double no_slope=Cst::PI);
void CalculateHick(double A[4][4], double& slope, double& Nx, double& Ny, double& Nz);
void CalculateFleming(double A[4][4], double& slope, double& Nx, double& Ny, double& Nz);
void CalculateHorn(double A[4][4], double& slope, double& Nx, double& Ny, double& Nz);
void CalculateCorripio(double A[4][4], double& slope, double& Nx, double& Ny, double& Nz);
void (DEMObject::*CalculateSlope)(double A[4][4], double& slope, double& Nx, double& Ny, double& Nz);
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);
double getCurvature(double A[4][4]);
double steepestGradient(double A[4][4]);
......
......@@ -40,13 +40,13 @@ double round(const double& x) {
namespace mio {
const int Date::daysLeapYear[12] = {31,29,31,30,31,30,31,31,30,31,30,31};
const int Date::daysNonLeapYear[12] = {31,28,31,30,31,30,31,31,30,31,30,31};
const double Date::DST_shift = 1.0; //in hours
const float Date::MJD_offset = 2400000.5; ///<offset between julian date and modified julian date
const float Date::Unix_offset = 2440587.5; ///<offset between julian date and Unix Epoch time
const float Date::Excel_offset = 2415018.5; ///<offset between julian date and Excel dates (note that excel invented some days...)
const float Date::Matlab_offset = 1721058.5; ///<offset between julian date and Matlab dates
const int Date::daysLeapYear[12] = {31,29,31,30,31,30,31,31,30,31,30,31};
const int Date::daysNonLeapYear[12] = {31,28,31,30,31,30,31,31,30,31,30,31};
const double Date::epsilon=1./(24.*3600.); ///< minimum difference between two dates. 1 second in units of days
//NOTE: For the comparison operators, we assume that dates are positive so we can bypass a call to abs()
......
......@@ -85,13 +85,13 @@ class Date {
CLOSEST ///< rounding toward closest
} RND;
static const int daysLeapYear[];
static const int daysNonLeapYear[];
static const double DST_shift;
static const float MJD_offset;
static const float Unix_offset;
static const float Excel_offset;
static const float Matlab_offset;
static const int daysLeapYear[];
static const int daysNonLeapYear[];
Date();
Date(const double& julian_in, const double& in_timezone, const bool& in_dst=false);
......
......@@ -98,8 +98,6 @@ class MeteoData : POPBase {
class MeteoData {
#endif
public:
static const size_t nrOfParameters; ///<holds the number of meteo parameters stored in MeteoData
/// \anchor meteoparam this enum provides indexed access to meteorological fields
enum Parameters {firstparam=0,
TA=firstparam, ///< Air temperature
......@@ -229,6 +227,8 @@ class MeteoData {
Date date; ///<Timestamp of the measurement
StationData meta; ///<The meta data of the measurement
static const size_t nrOfParameters; ///<holds the number of meteo parameters stored in MeteoData
private:
//static methods
static std::map<size_t, std::string> static_meteoparamname; ///<Associate a name with meteo parameters in Parameters
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment