diff --git a/src/eckit_readers/MvLocation.cc b/src/eckit_readers/MvLocation.cc index 9191cac0..04bc03b9 100644 --- a/src/eckit_readers/MvLocation.cc +++ b/src/eckit_readers/MvLocation.cc @@ -7,74 +7,142 @@ ***************************** LICENSE END *************************************/ -// MvLocation.cc, vk 940901... -// rev vk 011025 +#include +#include "MvLocation.h" -// classes: MvLocation, MvArea, MvXSectionLine +//=================================================== +// +// MvLocation +// +//=================================================== +const double MvLocation::cRadian = 180. / M_PI; +const double MvLocation::cDegree = 1 / MvLocation::cRadian; +const double MvLocation::cRadianToMetre = 1852 *180.0 * 60.0 / M_PI; +const double MvLocation::cMetreToRadian = 1./MvLocation::cRadianToMetre; -#include "MvLocation.h" -#include -//_____________________________________________________________________ set -void MvLocation ::set(double aLat, double aLong) { - fLatitude = aLat; - fLongitude = aLong; +MvLocation& MvLocation::operator=(const MvLocation& aLoc) +{ + set(aLoc.latitude(), aLoc.longitude()); + return *this; +} + +void MvLocation::set(double aLat, double aLong) +{ + lat_ = aLat; + lon_ = aLong; } -//______________________________________________________________ distanceInRadians -double MvLocation ::distanceInRadians(const MvLocation& anOtherLocation) const { - const double cDEG2RAD = M_PI / 180.0; +void MvLocation::ensureLongitudeBelow360() +{ + while (lon_ >= 360) { + lon_ -= 360; + } +} - double lat1 = latitude() * cDEG2RAD; - double lat2 = anOtherLocation.latitude() * cDEG2RAD; - double lon1 = longitude() * cDEG2RAD; - double lon2 = anOtherLocation.longitude() * cDEG2RAD; +// compute the cosine of the angular distance between the current and the other location +double MvLocation::cosOfDistance(double aLat, double aLon) const +{ + // short-circuit if the two points are the same, because floating-point + // imprecisions can cause a 'nan' on the next calculation + if ((lat_ == aLat) && (lon_ == aLon)) { + return 1.; + } - //-- from: http://williams.best.vwh.net/avform.htm (020815/vk) -- - double d = acos(sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(lon1 - lon2)); + double lat1 = degToRad(lat_); + double lat2 = degToRad(aLat); + double dlon = degToRad(lon_ - aLon); - return d; + //-- from: http://williams.best.vwh.net/avform.htm + return std::sin(lat1) * std::sin(lat2) + std::cos(lat1) * std::cos(lat2) * std::cos(dlon); } -//_____________________________________________________________ distanceInDegrees -double MvLocation ::distanceInDegrees(const MvLocation& anOtherLocation) const { - return distanceInRadians(anOtherLocation) * 180.0 / M_PI; +// compute the cosine of the angular distance between the current and the other location +double MvLocation::cosOfDistance(const MvLocation& other) const +{ + return cosOfDistance(other.latitude(), other.longitude()); } -//_____________________________________________________________ distanceInMeters -double MvLocation ::distanceInMeters(const MvLocation& anOtherLocation) const { - const double cR2NM = 180.0 * 60.0 / M_PI; //-- radians -> nautical miles - const double cNM2M = 1852; //-- nautical miles -> metres +// compute the angular distance in radians between the current and the other location +double MvLocation::distanceInRadians(const MvLocation& other) const +{ + return std::acos(cosOfDistance(other)); +} - double nm = distanceInRadians(anOtherLocation) * cR2NM; - return nm * cNM2M; - ; +// compute the angular distance in degrees between the current and the other location +double MvLocation::distanceInDegrees(const MvLocation& other) const +{ + return radToDeg(distanceInRadians(other)); } -//_____________________________________________________________ operator= -MvLocation& MvLocation ::operator=(const MvLocation& aLoc) { - set(aLoc.latitude(), aLoc.longitude()); - return *this; +// compute the distance in metres on the surface of Earth +// between the current and the other location +double MvLocation::distanceInMeters(const MvLocation& other) const +{ + return radiansToMetres(distanceInRadians(other)); } -//_____________________________________________________________ operator<< -ostream& operator<<(ostream& aStream, const MvLocation& aLocation) { +std::ostream& operator<<(std::ostream& aStream, const MvLocation& aLocation) +{ aStream << "(" << aLocation.latitude() << "," << aLocation.longitude() << ")"; + /*--- How to get a constant field width + constant nr of decimal digits!? + aStream << "("; + aStream.width( 6 ); aStream.fill( ' ' ); + aStream << aLocation.latitude() << ","; + aStream.width( 7 ); aStream.fill( ' ' ); + aStream << aLocation.longitude() << ")"; +---*/ return aStream; } +//=================================================== +// +// MvLocationHub +// +//=================================================== + +MvLocationHub& MvLocationHub::operator=(const MvLocationHub& aLoc) +{ + set(aLoc.latitude(), aLoc.longitude()); + cosLat_ = aLoc.cosLat_; + sinLat_ = aLoc.sinLat_; + return *this; +} + +double MvLocationHub::cosOfDistance(double aLat, double aLon) const +{ + if (sinLat_ < -100) { + sinLat_ = std::sin(degToRad(lat_)); + cosLat_ = std::cos(degToRad(lat_)); + } + + // short-circuit if the two points are the same, because floating-point + // imprecisions can cause a 'nan' on the next calculation + if (lat_ == aLat && lon_ == aLon) { + return 1.; + } + + double lat2 = degToRad(aLat); + double dLon = degToRad(lon_ - aLon); + + //-- from: http://williams.best.vwh.net/avform.htm + return sinLat_* std::sin(lat2) + cosLat_ * std::cos(lat2) * std::cos(dLon); +} + //____________________________________________________________________________ //============================================================================ MvArea //____________________________________________________________________________ -MvArea ::MvArea(void) { +MvArea ::MvArea() +{ MvLocation myLocation; set(myLocation, myLocation); } //____________________________________________________________________ -void MvArea ::set(const MvLocation& aLocation1, const MvLocation& aLocation2) { +void MvArea ::set(const MvLocation& aLocation1, const MvLocation& aLocation2) +{ double y1 = aLocation1.latitude(); double x1 = aLocation1.longitude(); double y2 = aLocation2.latitude(); @@ -94,22 +162,28 @@ void MvArea ::set(const MvLocation& aLocation1, const MvLocation& aLocation2) { } //____________________________________________________________________ -bool MvArea ::inside(const MvLocation& aPoint) const { - if (aPoint.latitude() >= fLowerLeft.latitude() && aPoint.latitude() <= fUpperRight.latitude() && - aPoint.longitude() >= fLowerLeft.longitude() && aPoint.longitude() <= fUpperRight.longitude()) +bool MvArea ::inside(const MvLocation& aPoint) const +{ + if (aPoint.latitude() >= fLowerLeft.latitude() && + aPoint.latitude() <= fUpperRight.latitude() && + aPoint.longitude() >= fLowerLeft.longitude() && + aPoint.longitude() <= fUpperRight.longitude()) return true; else return false; } //____________________________________________________________________ -MvArea& MvArea ::operator=(const MvArea& anArea) { +MvArea& +MvArea ::operator=(const MvArea& anArea) +{ set(anArea.lowerLeft(), anArea.upperRight()); return *this; } //____________________________________________________________________ -ostream& operator<<(ostream& aStream, const MvArea& anArea) { +std::ostream& operator<<(std::ostream& aStream, const MvArea& anArea) +{ aStream << anArea.lowerLeft() << "-" << anArea.upperRight(); return aStream; } @@ -119,7 +193,8 @@ ostream& operator<<(ostream& aStream, const MvArea& anArea) { //______________________________________________________________________ //_____________________________________________________________ WithinDelta -bool MvXSectionLine ::withinDelta(const MvLocation& aLocation) const { +bool MvXSectionLine ::withinDelta(const MvLocation& aLocation) const +{ //-- check that max delta has been set if (fMaxDeltaInMeters < 0) return false; @@ -132,7 +207,8 @@ bool MvXSectionLine ::withinDelta(const MvLocation& aLocation) const { return insideXLine(aLocation); } //_____________________________________________________________ InsideXLine -bool MvXSectionLine ::insideXLine(const MvLocation& aLocation) const { +bool MvXSectionLine ::insideXLine(const MvLocation& aLocation) const +{ MvLocation myLocation = nearestPointOnXLine(aLocation); MvArea myArea(fLocation1, fLocation2); return myArea.inside(myLocation); @@ -144,7 +220,9 @@ bool MvXSectionLine ::insideXLine(const MvLocation& aLocation) const { // WARNING: calculates distance from a line going through end // points of XSectionLine, NOT ONLY between points! //------------------------------------------------------------- -MvLocation MvXSectionLine ::nearestPointOnXLine(const MvLocation& aLocation) const { +MvLocation +MvXSectionLine ::nearestPointOnXLine(const MvLocation& aLocation) const +{ MvLocation myNearestPointOnXLine; double dy = fLocation1.latitude() - fLocation2.latitude(); @@ -177,19 +255,25 @@ MvLocation MvXSectionLine ::nearestPointOnXLine(const MvLocation& aLocation) con return myNearestPointOnXLine; } //_____________________________________________________________ DeltaInDegrees -double MvXSectionLine ::deltaInDegrees(const MvLocation& aLocation) const { +double +MvXSectionLine ::deltaInDegrees(const MvLocation& aLocation) const +{ return aLocation.distanceInDegrees(nearestPointOnXLine(aLocation)); } //_____________________________________________________________ DeltaInMeters // Q&D approximation for delta, uses DeltaInDegrees //------------------------------------------------------------- -double MvXSectionLine ::deltaInMeters(const MvLocation& aLocation) const { +double +MvXSectionLine ::deltaInMeters(const MvLocation& aLocation) const +{ //-- Q&D formula, out of my old used brains, unchecked!!!! (vk 940901) -- double degreeIntoMeters = 6370. * 1000. * 2. * M_PI / 360.; return deltaInDegrees(aLocation) * degreeIntoMeters; } //_____________________________________________________________ operator= -MvXSectionLine& MvXSectionLine ::operator=(const MvXSectionLine& anXLine) { +MvXSectionLine& +MvXSectionLine ::operator=(const MvXSectionLine& anXLine) +{ setLine(anXLine.startPoint(), anXLine.endPoint()); setMaxDelta(anXLine.maxDelta()); return *this; @@ -198,7 +282,9 @@ MvXSectionLine& MvXSectionLine ::operator=(const MvXSectionLine& anXLine) { // output format: "(lat1,long1)-(lat2,long2)/delta" //------------------------------------------------- -ostream& operator<<(ostream& aStream, const MvXSectionLine& anXLine) { - aStream << anXLine.startPoint() << "-" << anXLine.endPoint() << "/" << anXLine.fMaxDeltaInMeters; // << endl; +std::ostream& operator<<(std::ostream& aStream, const MvXSectionLine& anXLine) +{ + aStream << anXLine.startPoint() << "-" << anXLine.endPoint() + << "/" << anXLine.fMaxDeltaInMeters; // << std::endl; return aStream; } diff --git a/src/eckit_readers/MvLocation.h b/src/eckit_readers/MvLocation.h index a54bd1a5..81a5536b 100644 --- a/src/eckit_readers/MvLocation.h +++ b/src/eckit_readers/MvLocation.h @@ -7,85 +7,100 @@ ***************************** LICENSE END *************************************/ -// MvLocation.h, vk 940901... -// rev vk 950303 +#pragma once -#ifndef MvLocation_DEFINED_ -#define MvLocation_DEFINED_ - -#include "inc_iostream.h" +#include const double MISSING_LOC_VALUE = -99999.; +//========================================================= +// WARNING: we cannot use MvSci here because this code is +// also included in the Magics source!!! +//========================================================= -//_________________________________________________________________________ MvLocation //! Class for geographical locations /*! MvLocation is used to store latitude-longitude location values. * Class also provides methods to calculate the distance to another * geographical location */ -class MvLocation { +class MvLocation +{ //! Output operator for MvLocation object /*! The output is enclosed in parenthesis and latitude and longitude values - * are separated by a comma, for instance: - * 
-     *      MvLocation loc1(51.46,-1.33);
-     *      MvLocation loc2(60.45,25.0);
-     *      std::cout << "Locations are: " << loc1 << " and " << loc2 << std::endl;
-     *
- * would output the following line: - * 
-     *      Locations are: (51.46,-1.33) and (60.45,25)
-     *
- */ - friend ostream& operator<<(ostream& aStream, const MvLocation& aLocation); + * are separated by a comma, for instance: + * 
+ *      MvLocation loc1(51.46,-1.33);
+ *      MvLocation loc2(60.45,25.0);
+ *      std::cout << "Locations are: " << loc1 << " and " << loc2 << std::endl;
+ *
+ * would output the following line: + * 
+ *      Locations are: (51.46,-1.33) and (60.45,25)
+ *
+ */ + friend std::ostream& operator<<(std::ostream& aStream, const MvLocation& aLocation); public: - //! Constructor, location is assigned with missing values - MvLocation() { - fLatitude = MISSING_LOC_VALUE; - fLongitude = MISSING_LOC_VALUE; - } - - //! Constructor for latitude-longitude location + MvLocation() {} MvLocation(double aLat, double aLong) { set(aLat, aLong); } + MvLocation& operator=(const MvLocation& aLoc); - //! Sets new latitude-longitude location void set(double aLat, double aLong); - //! Checks that the stored location is a valid geographical point - /*! Latitude value must be in interval [-90,90] and longitude - * value in interval [-360,360], in degrees. - */ - bool ok() { return (fLatitude <= 90 && fLatitude >= -90 && fLongitude <= 360 && fLongitude >= -360); } - - //! Returns the latitude value - double latitude() const { return fLatitude; } - - //! Alias to method latitude() - double y() const { return fLatitude; } - - //! Returns the longitude value - double longitude() const { return fLongitude; } - - //! Alias to method longitude() - double x() const { return fLongitude; } - - //! Returns the distance (in radians) to the given point - double distanceInRadians(const MvLocation& anOtherLocation) const; - - //! Returns the distance (in degrees) to the given point - double distanceInDegrees(const MvLocation& anOtherLocation) const; - - //! Returns the distance (in metres) to the given point - double distanceInMeters(const MvLocation& anOtherLocation) const; + // Checks that the stored location is a valid geographical point + // Latitude value must be in interval [-90,90] and longitude + // value in interval [-360,360], in degrees. + bool ok() const { return (lat_ <= 90 && lat_ >= -90 && + lon_ <= 360 && lon_ >= -360); } + + // Subtracts 360 until the longitude is below 360 (e.g. 360 becomes 0, 370 becomes 10) + void ensureLongitudeBelow360(); + + double latitude() const { return lat_; } + double y() const { return lat_; } + double longitude() const { return lon_; } + double x() const { return lon_; } + + // Returns the cosine of the angular distance to the given point + double cosOfDistance(const MvLocation& other) const; + virtual double cosOfDistance(double aLat, double aLon) const; + + // Returns the distance (in radians) to the given point + double distanceInRadians(const MvLocation& other) const; + + // Returns the distance (in degrees) to the given point + double distanceInDegrees(const MvLocation& other) const; + + // Returns the distance (in metres) to the given point + double distanceInMeters(const MvLocation& other) const; + +protected: + // ideally McSci should be used here, but see warning above! + static double degToRad(double d) {return d * cDegree;} + static double radToDeg(double r) {return r * cRadian;} + static double metresToRadians(double m) {return m*cMetreToRadian;} + static double radiansToMetres(double r) {return r*cRadianToMetre;} + + static const double cRadian; + static const double cDegree; + static const double cMetreToRadian; + static const double cRadianToMetre; + double lat_{MISSING_LOC_VALUE}; + double lon_{MISSING_LOC_VALUE}; +}; - //! Assignment operator - MvLocation& operator=(const MvLocation& aLoc); +// A location from that the distance to other locations are computed. It +// caches some trigonometric values to make the computations faster. +class MvLocationHub : public MvLocation +{ +public: + using MvLocation::MvLocation; + MvLocationHub& operator=(const MvLocationHub& aLoc); + double cosOfDistance(double aLat, double aLon) const override; -private: - double fLatitude; - double fLongitude; +protected: + mutable double cosLat_{-1000.}; + mutable double sinLat_{-1000.}; }; //_________________________________________________________________________ MvArea @@ -93,12 +108,16 @@ class MvLocation { /*! This is another incarnation of MvGeoBox class, used mainly by * MvObsSetIterator. For other usage MvGeoBox is recommended over MvArea. */ -class MvArea { - friend ostream& operator<<(ostream& aStream, const MvArea& aArea); +class MvArea +{ + friend std::ostream& operator<<(std::ostream& aStream, const MvArea& aArea); public: MvArea(); - MvArea(const MvLocation& aLoc1, const MvLocation& aLoc2) { set(aLoc1, aLoc2); } + MvArea(const MvLocation& aLoc1, const MvLocation& aLoc2) + { + set(aLoc1, aLoc2); + } void set(const MvLocation& aLoc1, const MvLocation& aLoc2); bool inside(const MvLocation& aPoint) const; @@ -118,12 +137,14 @@ class MvArea { * distance (max delta) from the given line, i.e. close enough to be * considered to be within the line. */ -class MvXSectionLine { - friend ostream& operator<<(ostream& aStream, const MvXSectionLine& aXSectionLine); +class MvXSectionLine +{ + friend std::ostream& operator<<(std::ostream& aStream, const MvXSectionLine& aXSectionLine); public: //! Empty constructor creates a missing line - MvXSectionLine(void) { + MvXSectionLine(void) + { fLocation1.set(MISSING_LOC_VALUE, MISSING_LOC_VALUE); fLocation2.set(MISSING_LOC_VALUE, MISSING_LOC_VALUE); fMaxDeltaInMeters = -1; @@ -131,22 +152,25 @@ class MvXSectionLine { //! Constructor, only points defined, no max delta given /*! Use method setMaxDelta() to set the maximum allowed distance from the line - */ - MvXSectionLine(const MvLocation& aLoc1, const MvLocation& aLoc2) { + */ + MvXSectionLine(const MvLocation& aLoc1, const MvLocation& aLoc2) + { fLocation1 = aLoc1; fLocation2 = aLoc2; fMaxDeltaInMeters = -1; } //! Constructor, two points and max distance from the line - MvXSectionLine(const MvLocation& aLoc1, const MvLocation& aLoc2, double aDelta) { + MvXSectionLine(const MvLocation& aLoc1, const MvLocation& aLoc2, double aDelta) + { fLocation1 = aLoc1; fLocation2 = aLoc2; fMaxDeltaInMeters = aDelta; } //! Set a new line between points aLoc1 and aLoc2 - void setLine(const MvLocation& aLoc1, const MvLocation& aLoc2) { + void setLine(const MvLocation& aLoc1, const MvLocation& aLoc2) + { fLocation1 = aLoc1; fLocation2 = aLoc2; } @@ -188,5 +212,4 @@ class MvXSectionLine { double fMaxDeltaInMeters; }; -#endif // MvLocation_DEFINED_