WSL/SLF GitLab Repository

Commit 4b63f416 authored by Mathias Bavay's avatar Mathias Bavay
Browse files

A new group of classes has been created, under the meteolaws subdirectory....

A new group of classes has been created, under the meteolaws subdirectory. This contains various general meteorological laws, such as a sun radiation model, standard atmosphere, etc

The first implementation (not tested yet) of a non-linear least square fit algorithm has been implemented in the libfit1D. 

Several documentation issues have been fixed (obsolete code examples, etc) and classes grouped by modules. This should help the user to quickly locate the classes of interest. 
parent a7afd2a2
......@@ -10,6 +10,7 @@ PROJECT(meteoio)
ADD_SUBDIRECTORY(meteoio)
ADD_SUBDIRECTORY(meteoio/plugins)
ADD_SUBDIRECTORY(meteoio/plugins/gsn)
ADD_SUBDIRECTORY(meteoio/meteolaws)
# require proper c++
set(DEBUG_ARITHM ON CACHE BOOL "Use the debug flags -DDEBUG_ARITHM to catch arithmetic exceptions")
......
#include <stdio.h>
#include <stdlib.h>
#include <meteoio/MeteoIO.h>
int main(int /*argc*/, char** argv) {
mio::SunObject Sun(46.77181, 9.86820, 2192.); //Stillberg station
const double slope_azi=38., slope_elev=35.; //Stillberg station
const double TZ=1.;
const double TA = 273.15+11., RH = 0.5, mean_albedo = 0.5;
mio::Date d1;
d1.setTimeZone(TZ);
mio::IOUtils::convertString(d1,argv[1]); //get date and time
Sun.setDate(d1.getJulianDate(), d1.getTimeZone());
double iswr_ref;
mio::IOUtils::convertString(iswr_ref,argv[2]); //get measured global incoming radiation
Sun.calculateRadiation(TA, RH, mean_albedo);
std::cout << Sun;
//radiation in the beam
double B_toa, B_direct, B_diffuse;
Sun.getBeamRadiation(B_toa, B_direct, B_diffuse);
const double Md_beam=Sun.getSplitting(B_toa,iswr_ref);
//radiation on the horizontal
double H_toa, H_direct, H_diffuse;
Sun.getHorizontalRadiation(H_toa, H_direct, H_diffuse);
const double Md_horizontal=Sun.getSplitting(H_toa,iswr_ref);
//radiation on the slope
double S_toa, S_direct, S_diffuse;
Sun.getSlopeRadiation(slope_azi, slope_elev, S_toa, S_direct, S_diffuse);
const double Md_slope=Sun.getSplitting(S_toa,iswr_ref);
//outputing the results for beam, horizontal and slope
const unsigned int colw=10;
std::cout << std::setw(colw) << "Slope" << std::fixed << std::setw(colw) << std::setprecision(1) << S_toa;
std::cout << std::fixed << std::setw(colw) << std::setprecision(1) << S_direct;
std::cout << std::fixed << std::setw(colw) << std::setprecision(1) << S_diffuse;
std::cout << std::fixed << std::setw(colw) << std::setprecision(1) << S_direct+S_diffuse << "\n\n";
printf("Slope: azi=%g angle=%g ",slope_azi,slope_elev);
printf("Angle of incidence=%g\n",Sun.position.getAngleOfIncidence(slope_azi,slope_elev));
printf("Beam Splitting coefficient=%g\n",Md_beam);
printf("Horizontal Splitting coefficient=%g\n",Md_horizontal);
printf("Slope Splitting coefficient=%g\n",Md_slope);
return 0;
}
......@@ -32,6 +32,7 @@ namespace mio {
* @class Array
* @brief The template class Array is a 1D array (vector) able to hold any type of object as datatype
*
* @ingroup data_str
* @author Thomas Egger
* @date 2009-05-02
*/
......
......@@ -53,8 +53,10 @@ template <class T> class Array2DProxy {
/**
* @class Array2D
* @brief The template class Array2D is a 2D Array (Matrix) able to hold any type of object as datatype
* @brief The template class Array2D is a 2D Array (Matrix) able to hold any type of object as datatype.
* It relies on the Array2DProxy class to provide the [][] operator (slower than the (i,j) call).
*
* @ingroup data_str
* @author Thomas Egger
*/
template<class T> class Array2D {
......
......@@ -35,7 +35,7 @@ template <class T> class Array3DProxy2;
/**
* @class Array3DProxy
* @brief The template class Array3DProxy is a helper class for the template class Array3D
* with the purpose of adding the [][][] operator to Array3D
* with the purpose of adding the [][] operator to Array3D
*
* @author Thomas Egger
*/
......@@ -77,8 +77,9 @@ template <class T> class Array3DProxy2 {
/**
* @class Array3D
* @brief The template class Array3D is a 3D Array (Tensor) able to hold any type of object as datatype
*
* @brief The template class Array3D is a 3D Array (Tensor) able to hold any type of object as datatype.
* It relies on the Array3DProxy2 class to provide the [][][] operator (slower than the (i,j,k) call).
* @ingroup data_str
* @date 2009-07-19
* @author Thomas Egger
*/
......
......@@ -8,6 +8,8 @@ IF(PROJ4)
ADD_DEFINITIONS(-DPROJ4)
ENDIF(PROJ4)
INCLUDE("${PROJECT_SOURCE_DIR}/meteoio/meteolaws/CMakeLists.txt")
SET(meteoio_sources
IOManager.cc
ProcessingStack.cc
......@@ -35,6 +37,7 @@ SET(meteoio_sources
Grid3DObject.cc
IOInterface.cc
libinterpol1D.cc
libfit1D.cc
StationData.cc
Config.cc
DynamicLibrary.cc
......@@ -42,6 +45,7 @@ SET(meteoio_sources
IOUtils.cc
libinterpol2D.cc
MeteoData.cc
${meteolaws_sources}
)
if (DEST MATCHES "par")
......
......@@ -38,6 +38,8 @@ typedef void(Coords::*convfunc)(double, double, double&, double&) const;
* given object, as soon as a latitude/longitude can be computed/read, it will be used as a reference.
* This means that every subsequent change of projection system or read will be the conversion of this
* reference lat/lon position (until a new "set" is called). See Coords::setProj for the supported coordinate systems.
*
* @ingroup data_str
* @author Mathias Bavay
* @date 2009-01-20
*/
......
......@@ -36,6 +36,8 @@ class DEMObjectDummy {}; //HACK for POPC
* @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.
*
* @ingroup data_str
* @author Gaël Rosset - Mathias Bavay
* @date 2009-07-20
*/
......
......@@ -55,7 +55,8 @@ namespace mio {
* - truncated julian date, see Date::getTruncatedJulianDate
* - Unix date, see Date::getUnixDate
* - Excel date, see Date::getExcelDate
*
*
* @ingroup data_str
* @author Mathias Bavay
* @date 2010-04-15
*/
......
......@@ -35,7 +35,8 @@ namespace mio {
/**
* @class FilterAlgorithms
* @brief
* @brief A class to filter time series of data
*
* @author Thomas Egger
* @date 2009-11-03
*/
......
......@@ -32,6 +32,7 @@ namespace mio {
* @class Grid2DObject
* @brief A class to represent 2D Grids. Typical application as DEM or Landuse Model.
*
* @ingroup data_str
* @author Thomas Egger
* @date 2008-12-20
*/
......
......@@ -31,6 +31,7 @@ namespace mio {
* @class Grid3DObject
* @brief A class to represent 3D Grids. Typical application: wind field
*
* @ingroup data_str
* @author Thomas Egger
* @date 2009-07-20
*/
......
......@@ -84,6 +84,7 @@ namespace mio {
* a derived class is to be created that holds the specific implementation of the purely virtual member funtions.
* The IOHandler class is a wrapper class that is able to deal with all above implementations of the IOInterface abstract base class.
*
* @ingroup plugins
* @author Thomas Egger
* @date 2009-01-08
*/
......
......@@ -486,4 +486,18 @@ unsigned int IOUtils::seek(const Date& soughtdate, const std::vector<MeteoData>&
return IOUtils::npos;
}
std::string IOUtils::printFractionalDay(const double& fractional) {
const double hours=floor(fractional*24.);
const double minutes=floor((fractional*24.-hours)*60.);
const double seconds=fractional*24.*3600.-hours*3600.-minutes*60.;
std::stringstream tmp;
tmp << std::fixed << std::setfill('0') << std::setprecision(0);
tmp << std::setw(2) << hours << ":";
tmp << std::setw(2) << minutes << ":";
tmp << std::setw(2) << seconds;
return tmp.str();
}
} //namespace
......@@ -21,6 +21,7 @@
#include <meteoio/Coords.h>
#include <meteoio/Date.h>
#include <meteoio/IOExceptions.h>
#include <meteoio/meteolaws/Meteoconst.h>
#include <sstream>
#include <string>
......@@ -66,9 +67,8 @@ namespace IOUtils {
const short int snodata = -999;
const unsigned int npos = (unsigned int)-1; ///<npos is the out-of-range value
const double earth_radius = 6371e3; ///<Earth radius in meters
const double grid_epsilon = 5.; ///<What is an acceptable small distance on a grid, in meters
const double lon_epsilon = grid_epsilon / earth_radius; ///<in degrees. Small angle for longitudes, so sin(x)=x
const double lon_epsilon = grid_epsilon / Cst::earth_R0; ///<in degrees. Small angle for longitudes, so sin(x)=x
const double lat_epsilon = lon_epsilon/2.; ///<in degrees. Small angle for latitudes. Since for longitudes it is for 360deg, it has to be 180deg for latitudes
typedef enum NODATA_HANLDING {
......@@ -296,7 +296,15 @@ namespace IOUtils {
*/
void getTimeZoneParameters(const Config& cfg, double& tz_in, double& tz_out);
/**
* @brief Nicely format an hour given as fractional day into a human readable hour.
* @param fractional fractional day (ie: fractional part of a julian date)
* @return string containing a human readable time
*/
std::string printFractionalDay(const double& fractional);
} //end namespace IOUtils
} //end namespace mio
#endif
......@@ -16,6 +16,7 @@
along with MeteoIO. If not, see <http://www.gnu.org/licenses/>.
*/
#include <meteoio/InterpolationAlgorithms.h>
#include <meteoio/meteolaws/Atmosphere.h>
using namespace std;
......@@ -462,7 +463,7 @@ void RHAlgorithm::calculate(Grid2DObject& grid)
//Compute dew point temperatures at stations
for (unsigned int ii=0; ii<vecDataRH.size(); ii++){
vecTd[ii] = Interpol2D::RhtoDewPoint(vecDataRH[ii], vecDataTA[ii], 1);
vecTd[ii] = Atmosphere::RhtoDewPoint(vecDataRH[ii], vecDataTA[ii], 1);
}
//Krieging on Td
......@@ -476,7 +477,7 @@ void RHAlgorithm::calculate(Grid2DObject& grid)
//Recompute Rh from the interpolated td
for (unsigned int jj=0; jj<grid.nrows; jj++) {
for (unsigned int ii=0; ii<grid.ncols; ii++) {
grid.grid2D(ii,jj) = Interpol2D::DewPointtoRh(grid.grid2D(ii,jj), ta.grid2D(ii,jj), 1);
grid.grid2D(ii,jj) = Atmosphere::DewPointtoRh(grid.grid2D(ii,jj), ta.grid2D(ii,jj), 1);
}
}
}
......
......@@ -122,6 +122,8 @@ class Meteo2DInterpolator; // forward declaration, cyclic header include
/**
* @class InterpolationAlgorithm
* @brief A class to perform 2D spatial interpolations. For more, see \ref interpol2d
*
* @ingroup stats
* @author Thomas Egger
* @date 2010-04-01
*/
......
......@@ -19,6 +19,23 @@
#define __MAINPAGE_H__
namespace mio {
//groups
/*! \defgroup meteolaws Meteorological Laws
Documentation for meteorological laws and constants.
*/
/*! \defgroup plugins IO plugins
Documentation for available IO plugins. Please consider having a look at the \subpage plugins "Available plugins and usage" page and the \subpage dev_plugins "How to write a Plugin" page.
*/
/*! \defgroup stats Statistical calculations
Documentation for available statistical calculations. This is heavily used by the \subpage filters "Filters" as well as the \subpage resampling "1D interpolations" and \subpage interpol2d "2D interpolations".
*/
/*! \defgroup data_str Data structures
Documentation for available data structures.
*/
/**
* @mainpage Welcome to MeteoIO
* @section intro_sec Introduction
......@@ -296,19 +313,17 @@ namespace mio {
* Date d1;
* std::vector<mio::MeteoData> vecMeteo;
*
* mio::IOHandler *raw_io = NULL;
* mio::BufferedIOHandler *io = NULL;
* mio::IOManager *io = NULL;
*
* try {
* mio::Config cfg("io.ini");
* raw_io = new mio::IOHandler(cfg);
* io = new mio::BufferedIOHandler(*raw_io, cfg);
* io = new mio::IOManager(cfg);
* } catch (IOException& e){
* std::cout << "Problem with IOHandler creation, cause: " << e.what() << std::endl;
* }
*
* try {
* mio::convertString(d1,argv[1]);
* mio::IOUtils::convertString(d1,argv[1]);
* io->readMeteoData(d1, vecMeteo);
* } catch (IOException& e){
* std::cout << "Problem when reading data, cause: " << e.what() << std::endl;
......@@ -321,7 +336,6 @@ namespace mio {
* }
*
* delete io;
* delete raw_io;
*
* return 0;
* }
......@@ -334,17 +348,17 @@ namespace mio {
* int main(void) {
* const double dist_x=700, dist_y=1200;
* mio::DEMObject dem;
* mio::IOHandler *raw_io = NULL;
* mio::IOManager *io = NULL;
* mio::Config *cfg = NULL;
* int i,j;
*
* try {
* cfg = new mio::Config("io.ini");
* raw_io = new mio::IOHandler(cfg);
* io = new mio::IOManager(cfg);
* } catch (IOException& e){
* std::cout << "Problem with IOHandler creation, cause: " << e.what() << std::endl;
* }
* raw_io->readDEM(dem);
* io->readDEM(dem);
* mio::Coords point(*cfg);
* point.setLatLon(46.1592, 8.12993);
* dem.WGS84_to_grid(point, i,j);
......@@ -359,18 +373,31 @@ namespace mio {
* }
* \endcode
*
* The next example shows how to compute and output spatial interpolations using previously read meteorological data and DEM.
* The next example shows how to compute and output spatial interpolations.
* \code
* void spatial_interpolations(mio::IOHandler& io, mio::DEMObject& dem, std::vector<mio::MeteoData>& vecMeteo);
* {
* mio::Grid2DObject ta;
* mio::Meteo2DInterpolator mi(dem, vecMeteo);
*
* mi.interpolate(mio::MeteoData::TA, ta);
*
* std::cout << "Writing the Grids to files" << std::endl;
* io->write2DGrid(ta, "ta_2d.asc");
* #include "MeteoIO.h"
*
* int main(void) {
* mio::Date d1;
*
* //initializing the io handlers according to the config file
* mio::Config cfg("io.ini");
* mio::IOManager io(cfg);
*
* //reading the dem (necessary for several spatial interpolations algoritms)
* mio::DEMObject dem;
* io.readDEM(dem);
*
* //we assume that the time given on the command line is in TZ=+1
* d1.setTimeZone(1.);
* mio::IOUtils::convertString(d1,argv[1]);
*
* //performing spatial interpolations
* mio::Grid2DObject ta_grid;
* io.interpolate(d1, dem, MeteoData::TA, ta_grid);
* io.write2DGrid(param,"ta.asc");
*
* return 0;
* }
* \endcode
*/
......
......@@ -335,6 +335,34 @@ double Matrix::scalar() const {
return operator()(1,1);
}
double Matrix::dot(const Matrix& A, const Matrix& B) {
unsigned int Acols, Arows, Bcols, Brows;
A.size(Arows, Acols);
B.size(Brows, Bcols);
if(Acols!=1 || Bcols!=1) {
std::stringstream tmp;
tmp << "Trying to get dot product of non vector matrix ";
tmp << "(" << Arows << "," << Acols << ") · ";
tmp << "(" << Brows << "," << Bcols << ") · ";
throw IOException(tmp.str(), AT);
}
if(Arows!=Brows) {
std::stringstream tmp;
tmp << "Trying to get dot product of incompatible matrix ";
tmp << "(" << Arows << "," << Acols << ") · ";
tmp << "(" << Brows << "," << Bcols << ") · ";
throw IOException(tmp.str(), AT);
}
double sum=0.;
for(unsigned int i=1; i<=Arows; i++) {
sum += A(i,1)*B(i,1);
}
return sum;
}
Matrix Matrix::T(const Matrix& m) {
return m.getT();
}
......
......@@ -38,6 +38,7 @@ namespace mio {
* It might be later possible to chose between using the embedded implementation or to act as a
* front end to BLAS for those who have it installed on their system.
*
* @ingroup data_str
* @author Mathias Bavay
*/
class Matrix {
......@@ -111,6 +112,12 @@ class Matrix {
double scalar() const;
static double scalar(const Matrix& m);
/**
* @brief Dot product.
* @return scalar value
*/
static double dot(const Matrix& A, const Matrix& B);
/**
* @brief matrix transpose.
* @return transposed matrix
......
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