WSL/SLF GitLab Repository

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

Following an issue with an wrongly chosen coordinate system, a few...

Following an issue with an wrongly chosen coordinate system, a few improvements have been developed: the documentation has been improved, UTM zones can be provided as {num}{letter} or the opposite, and a bug has been fixed in NetCDFIO (in some cases, the opened file would not be closed).
parent 4583c590
......@@ -377,7 +377,7 @@ void Meteo2DInterpolator::initVirtualStations(const bool& adjust_coordinates)
if (!dem.gridify(curr_point)) {
ostringstream ss;
ss << "Virtual station \"" << vecStation[ii] << "\" is not contained is provided DEM";
ss << "Virtual station \"" << vecStation[ii] << "\" is not contained is provided DEM " << dem.toString(DEMObject::SHORT);
throw NoDataException(ss.str(), AT);
}
......
......@@ -40,12 +40,13 @@ namespace mio {
* geolocalized data, the desired coordinate system must also be specified for the outputs (in the output section).
* This is done through the use of the COORDIN and COORDPARAM keys (see the documentation for each plugin).
*
* \anchor Coordinate_types
* There are two ways of supporting a given coordinate system: through the use of an adhoc implementation
* (that becomes part of MeteoIO) or through the use of an external library, Proj4 [ref: http://trac.osgeo.org/proj/].
* The current internal implementations are the following (given by their keyword):
* - CH1903 or CH1903+ for coordinates in the Swiss Grid [ref: http://geomatics.ladetto.ch/ch1903_wgs84_de.pdf]
* - UTM for UTM coordinates (the zone must be specified in the parameters, for example 31T) [ref: http://www.oc.nps.edu/oc2902w/maps/utmups.pdf]
* - UPS for Universal Polar Stereographic coordinates (the zone, either N or S, must be specified in the parameters). [ref: J. Hager, J. Behensky, B. Drew, <i>THE UNIVERSAL GRIDS: Universal Transverse Mercator (UTM) and Universal Polar Stereographic (UPS)</i>, 1989, Defense Mapping Agency, DMATM 8358.2]
* - <a href="https://en.wikipedia.org/wiki/Swiss_coordinate_system">CH1903 or CH1903+</a> for coordinates in the Swiss Grid [ref: http://geomatics.ladetto.ch/ch1903_wgs84_de.pdf]
* - <a href="https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system">UTM</a> for UTM coordinates (the zone must be specified in the parameters, for example 31T) [ref: http://www.oc.nps.edu/oc2902w/maps/utmups.pdf]
* - <a href="https://en.wikipedia.org/wiki/Universal_polar_stereographic_coordinate_system">UPS</a> for Universal Polar Stereographic coordinates (the zone, either N or S, must be specified in the parameters). [ref: J. Hager, J. Behensky, B. Drew, <i>THE UNIVERSAL GRIDS: Universal Transverse Mercator (UTM) and Universal Polar Stereographic (UPS)</i>, 1989, Defense Mapping Agency, DMATM 8358.2]
* - LOCAL for local coordinate system (using the horizontal and vertical distance from a reference point, see Coords::geo_distances for the available choice of distance algorithms)
*
* Such an example of use is the following:
......@@ -54,7 +55,7 @@ namespace mio {
* COORDPARAM = 31T
* @endcode
*
* On the other hand, when using the Proj4 library for handling the coordinate conversion, the EPSG codes of
* On the other hand, when using the <a href="https://en.wikipedia.org/wiki/PROJ.4">Proj4</a> library for handling the coordinate conversion, the EPSG codes of
* the chosen projection must be specified (such codes can be found at http://spatialreference.org/ref/epsg/?page=1)
* as illustrated below (21781 is the EPSG code for the CH1903 coordinate system. Such a code is 32767 at the maximum):
* @code
......@@ -596,16 +597,9 @@ void Coords::setAltitude(const double in_altitude, const bool in_update) {
/**
* @brief Set projection to use
* This projection will be used for converting between lat/lon and East/North
* This projection will be used for converting between lat/lon and East/North (see the \ref Coordinate_types "supported projections")
* @param[in] in_coordinatesystem string identifying the coordinate system to use
* @param[in] in_parameters string giving some additional parameters for the projection (optional)
*
* \anchor Coordinate_types
* The coordinate system can be any of the following:
* - CH1903 for coordinates in the Swiss Grid [ref: http://geomatics.ladetto.ch/ch1903_wgs84_de.pdf]
* - UTM for UTM coordinates (the zone must be specified in the parameters, for example 31T) [ref: http://www.oc.nps.edu/oc2902w/maps/utmups.pdf]
* - PROJ4 for coordinate conversion relying on the Proj4 library [ref: http://trac.osgeo.org/proj/]
* - LOCAL for local coordinate system (using the horizontal and vertical distance from a reference point, see Coords::geo_distances for the available choice of distance algorithms)
*/
void Coords::setProj(const std::string& in_coordinatesystem, const std::string& in_parameters)
{
......
......@@ -710,7 +710,9 @@ void CoordsAlgorithms::UPS_to_WGS84(const double& east_in, const double& north_i
void CoordsAlgorithms::parseUTMZone(const std::string& zone_info, char& zoneLetter, short int& zoneNumber)
{ //helper method: parse a UTM zone specification string into letter and number
if ((sscanf(zone_info.c_str(), "%hd%c", &zoneNumber, &zoneLetter) < 2) &&
(sscanf(zone_info.c_str(), "%hd %c)", &zoneNumber, &zoneLetter) < 2)) {
(sscanf(zone_info.c_str(), "%hd %c)", &zoneNumber, &zoneLetter) < 2) &&
(sscanf(zone_info.c_str(), "%c%hd", &zoneLetter, &zoneNumber) < 2) &&
(sscanf(zone_info.c_str(), "%c %hd)", &zoneLetter, &zoneNumber) < 2)) {
throw InvalidFormatException("Can not parse given UTM zone: "+zone_info,AT);
}
if (zoneNumber<1 || zoneNumber>60) {
......
......@@ -916,17 +916,35 @@ double DEMObject::safeGet(const int i, const int j)
return grid2D((unsigned)i, (unsigned)j);
}
const std::string DEMObject::toString() const {
const std::string DEMObject::toString(const FORMATS& type) const
{
std::ostringstream os;
os << "<DEMObject>\n";
os << llcorner.toString();
os << grid2D.getNx() << " x " << grid2D.getNy() << " @ " << cellsize << "m\t[ " << min_altitude << " - " << max_altitude << " ]\n";
//os << grid2D.toString();
os << "Slope: " << slope.getNx() << " x " << slope.getNy() << "\t[ " << min_slope << " - " << max_slope << " ]\n";
os << "Azi: " << azi.getNx() << " x " << azi.getNy() << "\n";
os << "Curvature: " << curvature.getNx() << " x " << curvature.getNy()<< "\t[ " << min_curvature << " - " << max_curvature << " ]\n";
os << "Nx: " << Nx.getNx() << " x " << Nx.getNy() << " , Ny: " << Ny.getNx() << " x " << Ny.getNy() << " , Nz: " << Nz.getNx() << " x " << Nz.getNy() << "\n";
os << "</DEMObject>\n";
switch(type) {
case(FULL):
{
os << "<DEMObject>\n";
os << llcorner.toString();
os << grid2D.getNx() << " x " << grid2D.getNy() << " @ " << cellsize << "m\t[ " << min_altitude << " - " << max_altitude << " ]\n";
//os << grid2D.toString();
os << "Slope: " << slope.getNx() << " x " << slope.getNy() << "\t[ " << min_slope << " - " << max_slope << " ]\n";
os << "Azi: " << azi.getNx() << " x " << azi.getNy() << "\n";
os << "Curvature: " << curvature.getNx() << " x " << curvature.getNy()<< "\t[ " << min_curvature << " - " << max_curvature << " ]\n";
os << "Nx: " << Nx.getNx() << " x " << Nx.getNy() << " , Ny: " << Ny.getNx() << " x " << Ny.getNy() << " , Nz: " << Nz.getNx() << " x " << Nz.getNy() << "\n";
os << "</DEMObject>\n";
break;
}
case(SHORT):
{
Coords urcorner(llcorner);
urcorner.moveByXY(static_cast<double>(grid2D.getNx())*cellsize, static_cast<double>(grid2D.getNy())*cellsize);
os << llcorner.toString(Coords::LATLON) << " / " << urcorner.toString(Coords::LATLON);
break;
}
default:
throw InvalidFormatException("Unsupported DEM information format", AT);
}
return os.str();
}
......
......@@ -63,6 +63,12 @@ class DEMObject : public Grid2DObject {
CURVATURE=4 ///< update the curvatures
} update_type;
///Keywords for selecting the toString formats
typedef enum {
FULL, ///< Provide all the usually necessary information
SHORT ///< Simplified, lat/lon only
} FORMATS;
DEMObject(const slope_type& i_algorithm=DFLT);
DEMObject(const size_t& ncols_in, const size_t& nrows_in,
......@@ -121,7 +127,7 @@ class DEMObject : public Grid2DObject {
bool operator==(const DEMObject& in) const; ///<Operator that tests for equality
bool operator!=(const DEMObject& in) const; ///<Operator that tests for inequality
const std::string toString() const;
const std::string toString(const FORMATS& type = FULL) const;
friend std::iostream& operator<<(std::iostream& os, const DEMObject& dem);
friend std::iostream& operator>>(std::iostream& is, DEMObject& dem);
......
......@@ -546,7 +546,10 @@ bool NetCDFIO::read2DGrid_internal(Grid2DObject& grid_out, const std::string& fi
int ncid, varid;
ncpp::open_file(filename, NC_NOWRITE, ncid);
if (!ncpp::check_variable(ncid, varname)) return false;
if (!ncpp::check_variable(ncid, varname)) {
ncpp::close_file(filename, ncid);
return false;
}
ncpp::get_variable(ncid, varname, varid);
std::vector<int> dimid, dim_varid;
......
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