36 #define M_PI 3.14159265358979323846
39 #define DEG2RAD(x) ((x)*M_PI/180)
45 mEllipsoidalMode =
false;
60 delete mCoordTransform;
66 if (
this == & origDA )
78 mEllipsoidalMode = origDA.mEllipsoidalMode;
79 mEllipsoid = origDA.mEllipsoid;
80 mSemiMajor = origDA.mSemiMajor;
81 mSemiMinor = origDA.mSemiMinor;
82 mInvFlattening = origDA.mInvFlattening;
91 mEllipsoidalMode = flag;
121 sqlite3_stmt *myPreparedStatement;
138 bool semiMajorOk, semiMinorOk;
139 double semiMajor = paramList[1].toDouble( & semiMajorOk );
140 double semiMinor = paramList[2].toDouble( & semiMinorOk );
141 if ( semiMajorOk && semiMinorOk )
163 QString mySql =
"select radius, parameter2 from tbl_ellipsoid where acronym='" + ellipsoid +
"'";
164 myResult = sqlite3_prepare( myDatabase, mySql.
toUtf8(), mySql.
toUtf8().
length(), &myPreparedStatement, &myTail );
166 if ( myResult == SQLITE_OK )
168 if ( sqlite3_step( myPreparedStatement ) == SQLITE_ROW )
170 radius =
QString((
char * )sqlite3_column_text( myPreparedStatement, 0 ) );
171 parameter2 =
QString((
char * )sqlite3_column_text( myPreparedStatement, 1 ) );
175 sqlite3_finalize( myPreparedStatement );
176 sqlite3_close( myDatabase );
181 QgsDebugMsg(
QString(
"setEllipsoid: no row in tbl_ellipsoid for acronym '%1'" ).arg( ellipsoid ) );
186 if ( radius.
left( 2 ) ==
"a=" )
190 QgsDebugMsg(
QString(
"setEllipsoid: wrong format of radius field: '%1'" ).arg( radius ) );
197 if ( parameter2.
left( 2 ) ==
"b=" )
200 mInvFlattening = mSemiMajor / ( mSemiMajor - mSemiMinor );
202 else if ( parameter2.
left( 3 ) ==
"rf=" )
205 mSemiMinor = mSemiMajor - ( mSemiMajor / mInvFlattening );
209 QgsDebugMsg(
QString(
"setEllipsoid: wrong format of parameter2 field: '%1'" ).arg( parameter2 ) );
213 QgsDebugMsg(
QString(
"setEllipsoid: a=%1, b=%2, 1/f=%3" ).arg( mSemiMajor ).arg( mSemiMinor ).arg( mInvFlattening ) );
217 QString proj4 =
"+proj=longlat +ellps=" + ellipsoid +
" +no_defs";
224 if ( destCRS.
srsid() == 0 )
227 .
arg(
QObject::tr(
"Generated CRS",
"A CRS automatically generated from layer info get this prefix for description" ) )
249 mEllipsoid =
QString(
"PARAMETER:%1:%2" ).
arg( semiMajor ).
arg( semiMinor );
250 mSemiMajor = semiMajor;
251 mSemiMinor = semiMinor;
252 mInvFlattening = mSemiMajor / ( mSemiMajor - mSemiMinor );
264 const unsigned char* wkb = geometry->
asWkb();
273 double res, resTotal = 0;
277 bool hasZptr =
false;
294 for ( i = 0; i < count; i++ )
315 for ( i = 0; i < count; i++ )
339 const unsigned char* wkb = geometry->
asWkb();
347 double res = 0.0, resTotal = 0.0;
351 bool hasZptr =
false;
374 for ( i = 0; i < count; i++ )
405 for (
int i = 0; i < nPoints; ++i )
411 wkbPtr +=
sizeof( double );
423 if ( points.
size() < 2 )
431 if ( mEllipsoidalMode && ( mEllipsoid !=
GEO_NONE ) )
432 p1 = mCoordTransform->
transform( points[0] );
438 if ( mEllipsoidalMode && ( mEllipsoid !=
GEO_NONE ) )
472 if ( mEllipsoidalMode && ( mEllipsoid !=
GEO_NONE ) )
530 for (
int idx = 0; idx < numRings; idx++ )
537 for (
int jdx = 0; jdx < nPoints; jdx++ )
543 wkbPtr +=
sizeof( double );
548 if ( mEllipsoidalMode && ( mEllipsoid !=
GEO_NONE ) )
555 if ( points.
size() > 2 )
598 if ( mEllipsoidalMode && ( mEllipsoid !=
GEO_NONE ) )
626 if ( mEllipsoidalMode && ( mEllipsoid !=
GEO_NONE ) )
634 double dx = p2.
x() - p1.
x();
635 double dy = p2.
y() - p1.
y();
636 bearing = atan2( dx, dy );
648 double* course1,
double* course2 )
const
650 if ( p1.
x() == p2.
x() && p1.
y() == p2.
y() )
654 double a = mSemiMajor;
655 double b = mSemiMinor;
656 double f = 1 / mInvFlattening;
661 double L = p2_lon - p1_lon;
662 double U1 = atan(( 1 - f ) * tan( p1_lat ) );
663 double U2 = atan(( 1 - f ) * tan( p2_lat ) );
664 double sinU1 = sin( U1 ), cosU1 = cos( U1 );
665 double sinU2 = sin( U2 ), cosU2 = cos( U2 );
667 double lambdaP = 2 *
M_PI;
669 double sinLambda = 0;
670 double cosLambda = 0;
675 double cosSqAlpha = 0;
676 double cos2SigmaM = 0;
682 while ( qAbs( lambda - lambdaP ) > 1e-12 && --iterLimit > 0 )
684 sinLambda = sin( lambda );
685 cosLambda = cos( lambda );
686 tu1 = ( cosU2 * sinLambda );
687 tu2 = ( cosU1 * sinU2 - sinU1 * cosU2 * cosLambda );
688 sinSigma = sqrt( tu1 * tu1 + tu2 * tu2 );
689 cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
690 sigma = atan2( sinSigma, cosSigma );
691 alpha = asin( cosU1 * cosU2 * sinLambda / sinSigma );
692 cosSqAlpha = cos( alpha ) * cos( alpha );
693 cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
694 C = f / 16 * cosSqAlpha * ( 4 + f * ( 4 - 3 * cosSqAlpha ) );
696 lambda = L + ( 1 - C ) * f * sin( alpha ) *
697 ( sigma + C * sinSigma * ( cos2SigmaM + C * cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) ) );
700 if ( iterLimit == 0 )
703 double uSq = cosSqAlpha * ( a * a - b * b ) / ( b * b );
704 double A = 1 + uSq / 16384 * ( 4096 + uSq * ( -768 + uSq * ( 320 - 175 * uSq ) ) );
705 double B = uSq / 1024 * ( 256 + uSq * ( -128 + uSq * ( 74 - 47 * uSq ) ) );
706 double deltaSigma = B * sinSigma * ( cos2SigmaM + B / 4 * ( cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) -
707 B / 6 * cos2SigmaM * ( -3 + 4 * sinSigma * sinSigma ) * ( -3 + 4 * cos2SigmaM * cos2SigmaM ) ) );
708 double s = b * A * ( sigma - deltaSigma );
712 *course1 = atan2( tu1, tu2 );
717 *course2 = atan2( cosU1 * sinLambda, -sinU1 * cosU2 + cosU1 * sinU2 * cosLambda ) +
M_PI;
725 return sqrt(( p2.
x() - p1.
x() ) * ( p2.
x() - p1.
x() ) + ( p2.
y() - p1.
y() ) * ( p2.
y() - p1.
y() ) );
730 if ( points.
size() < 2 )
743 if ( mEllipsoidalMode && ( mEllipsoid !=
GEO_NONE ) )
772 double QgsDistanceArea::getQ(
double x )
const
779 return sinx *( 1 + sinx2 *( m_QA + sinx2 *( m_QB + sinx2 * m_QC ) ) );
783 double QgsDistanceArea::getQbar(
double x )
const
790 return cosx *( m_QbarA + cosx2 *( m_QbarB + cosx2 *( m_QbarC + cosx2 * m_QbarD ) ) );
802 double a2 = ( mSemiMajor * mSemiMajor );
803 double e2 = 1 - ( a2 / ( mSemiMinor * mSemiMinor ) );
811 m_AE = a2 * ( 1 - e2 );
813 m_QA = ( 2.0 / 3.0 ) * e2;
814 m_QB = ( 3.0 / 5.0 ) * e4;
815 m_QC = ( 4.0 / 7.0 ) * e6;
817 m_QbarA = -1.0 - ( 2.0 / 3.0 ) * e2 - ( 3.0 / 5.0 ) * e4 - ( 4.0 / 7.0 ) * e6;
818 m_QbarB = ( 2.0 / 9.0 ) * e2 + ( 2.0 / 5.0 ) * e4 + ( 4.0 / 7.0 ) * e6;
819 m_QbarC = - ( 3.0 / 25.0 ) * e4 - ( 12.0 / 35.0 ) * e6;
820 m_QbarD = ( 4.0 / 49.0 ) * e6;
822 m_Qp = getQ( M_PI / 2 );
823 m_E = 4 * M_PI * m_Qp * m_AE;
831 double x1, y1, x2, y2, dx, dy;
836 if (( ! mEllipsoidalMode ) || ( mEllipsoid ==
GEO_NONE ) )
840 int n = points.
size();
841 x2 =
DEG2RAD( points[n-1].x() );
842 y2 =
DEG2RAD( points[n-1].y() );
843 Qbar2 = getQbar( y2 );
847 for (
int i = 0; i < n; i++ )
855 Qbar2 = getQbar( y2 );
858 while ( x1 - x2 >
M_PI )
861 while ( x2 - x1 >
M_PI )
865 area += dx * ( m_Qp - getQ( y2 ) );
867 if (( dy = y2 - y1 ) != 0.0 )
868 area += dx * getQ( y2 ) - ( dx / dy ) * ( Qbar2 - Qbar1 );
870 if (( area *= m_AE ) < 0.0 )
880 if ( area > m_E / 2 )
892 size = points.
size();
895 for ( i = 0; i < size; i++ )
900 area = area + points[i].x() * points[( i+1 ) % size].y() - points[( i+1 ) % size].x() * points[i].y();
920 else if ( qAbs( value ) > 1000000.0 )
923 value = value / 1000000.0;
925 else if ( qAbs( value ) > 10000.0 )
928 value = value / 10000.0;
937 if ( keepBaseUnit || qAbs( value ) == 0.0 )
941 else if ( qAbs( value ) > 1000.0 )
944 value = value / 1000;
946 else if ( qAbs( value ) < 0.01 )
949 value = value * 1000;
951 else if ( qAbs( value ) < 0.1 )
965 if ( keepBaseUnit || qAbs( value ) <= 0.5*43560.0 )
970 else if ( qAbs( value ) <= 0.5*5280.0*5280.0 )
980 value /= 5280.0 * 5280.0;
985 if ( qAbs( value ) <= 528.0 || keepBaseUnit )
987 if ( qAbs( value ) == 1.0 )
1020 if ( qAbs( value ) == 1.0 )
1030 QgsDebugMsg(
QString(
"Error: not picked up map units - actual value = %1" ).arg( u ) );
1047 QgsDebugMsg(
"We're measuring on an ellipsoid or using projections, the system is returning meters" );
1049 else if ( mEllipsoidalMode && mEllipsoid ==
GEO_NONE )
1053 QgsDebugMsg(
"We're measuing on planimetric distance/area on given CRS, measured value is in CRS units" );
1059 factorUnits *= factorUnits;
1062 measure *= factorUnits;
1064 measureUnits = displayUnits;
double bearing(const QgsPoint &p1, const QgsPoint &p2) const
compute bearing - in radians
void convertMeasurement(double &measure, QGis::UnitType &measureUnits, QGis::UnitType displayUnits, bool isArea) const
Helper for conversion between physical units.
QString toString(qlonglong i) const
~QgsDistanceArea()
Destructor.
bool saveAsUserCRS(QString name)
Copied from QgsCustomProjectionDialog /// Please refactor into SQL handler !!! ///.
QStringList split(const QString &sep, SplitBehavior behavior, Qt::CaseSensitivity cs) const
void setSourceCrs(long srsid)
sets source spatial reference system (by QGIS CRS)
void computeAreaInit()
precalculates some values (must be called always when changing ellipsoid)
A geometry is the spatial representation of a feature.
WkbType
Used for symbology operations.
bool setEllipsoid(const QString &ellipsoid)
sets ellipsoid by its acronym
double toDouble(bool *ok) const
QString tr(const char *sourceText, const char *disambiguation, int n)
double measurePolygon(const QList< QgsPoint > &points) const
measures polygon area
QString trUtf8(const char *sourceText, const char *disambiguation, int n)
const QString GEO_NONE
Constant that holds the string representation for "No ellips/No CRS".
static void logMessage(QString message, QString tag=QString::null, MessageLevel level=WARNING)
add a message to the instance (and create it if necessary)
bool createFromOgcWmsCrs(QString theCrs)
Set up this CRS from the given OGC CRS.
double measure(const QgsGeometry *geometry) const
general measurement (line distance or polygon area)
QString number(int n, int base)
bool createFromSrsId(const long theSrsId)
double computeDistance(const QList< QgsPoint > &points) const
calculate distance with given coordinates (does not do a transform anymore)
void append(const T &value)
#define QgsDebugMsgLevel(str, level)
double computePolygonArea(const QList< QgsPoint > &points) const
calculates area of polygon on ellipsoid algorithm has been taken from GRASS: gis/area_poly1.c
double computeDistanceBearing(const QgsPoint &p1, const QgsPoint &p2, double *course1=NULL, double *course2=NULL) const
calculates distance from two points on ellipsoid based on inverse Vincenty's formulae ...
double computePolygonFlatArea(const QList< QgsPoint > &points) const
bool startsWith(const QString &s, Qt::CaseSensitivity cs) const
const long GEOCRS_ID
Magic number for a geographic coord sys in QGIS srs.db tbl_srs.srs_id.
double measureLine(const QList< QgsPoint > &points) const
measures line
QString toString() const
String representation of the point (x,y)
double measurePerimeter(const QgsGeometry *geometry) const
measures perimeter of polygon
A class to represent a point.
QgsDistanceArea & operator=(const QgsDistanceArea &origDA)
Assignment operator.
static QString textUnit(double value, int decimals, QGis::UnitType u, bool isArea, bool keepBaseUnit=false)
General purpose distance and area calculator.
QString mid(int position, int n) const
const QString & ellipsoid() const
returns ellipsoid's acronym
Class for storing a coordinate reference system (CRS)
static const QString srsDbFilePath()
Returns the path to the srs.db file.
static double fromUnitToUnitFactor(QGis::UnitType fromUnit, QGis::UnitType toUnit)
Returns the conversion factor between the specified units.
UnitType
Map units that qgis supports.
QString left(int n) const
Custom exception class for Coordinate Reference System related exceptions.
static QString toLiteral(QGis::UnitType unit)
Provides the canonical name of the type value.
QString arg(qlonglong a, int fieldWidth, int base, const QChar &fillChar) const
const unsigned char * asWkb() const
Returns the buffer containing this geometry in WKB format.
QgsDistanceArea()
Constructor.
double computeDistanceFlat(const QgsPoint &p1, const QgsPoint &p2) const
uses flat / planimetric / Euclidean distance
void setEllipsoidalMode(bool flag)
sets whether coordinates must be projected to ellipsoid before measuring
bool createFromProj4(const QString &theProjString)
QString toProj4() const
Get the Proj Proj4 string representation of this srs.
void setSourceAuthId(QString authid)
sets source spatial reference system by authid
QByteArray toUtf8() const
QGis::UnitType mapUnits() const