41 #define M_PI 3.14159265358979323846
44 #define DEG2RAD(x) ((x)*M_PI/180)
50 mEllipsoidalMode =
false;
59 : mCoordTransform( 0 )
66 delete mCoordTransform;
72 if (
this == & origDA )
84 mEllipsoidalMode = origDA.mEllipsoidalMode;
85 mEllipsoid = origDA.mEllipsoid;
86 mSemiMajor = origDA.mSemiMajor;
87 mSemiMinor = origDA.mSemiMinor;
88 mInvFlattening = origDA.mInvFlattening;
92 delete mCoordTransform;
98 mEllipsoidalMode = flag;
128 sqlite3_stmt *myPreparedStatement;
145 bool semiMajorOk, semiMinorOk;
146 double semiMajor = paramList[1].toDouble( & semiMajorOk );
147 double semiMinor = paramList[2].toDouble( & semiMinorOk );
148 if ( semiMajorOk && semiMinorOk )
170 QString mySql =
"select radius, parameter2 from tbl_ellipsoid where acronym='" + ellipsoid +
"'";
171 myResult = sqlite3_prepare( myDatabase, mySql.
toUtf8(), mySql.
toUtf8().
length(), &myPreparedStatement, &myTail );
173 if ( myResult == SQLITE_OK )
175 if ( sqlite3_step( myPreparedStatement ) == SQLITE_ROW )
177 radius =
QString((
char * )sqlite3_column_text( myPreparedStatement, 0 ) );
178 parameter2 =
QString((
char * )sqlite3_column_text( myPreparedStatement, 1 ) );
182 sqlite3_finalize( myPreparedStatement );
183 sqlite3_close( myDatabase );
188 QgsDebugMsg(
QString(
"setEllipsoid: no row in tbl_ellipsoid for acronym '%1'" ).arg( ellipsoid ) );
193 if ( radius.
left( 2 ) ==
"a=" )
197 QgsDebugMsg(
QString(
"setEllipsoid: wrong format of radius field: '%1'" ).arg( radius ) );
204 if ( parameter2.
left( 2 ) ==
"b=" )
207 mInvFlattening = mSemiMajor / ( mSemiMajor - mSemiMinor );
209 else if ( parameter2.
left( 3 ) ==
"rf=" )
212 mSemiMinor = mSemiMajor - ( mSemiMajor / mInvFlattening );
216 QgsDebugMsg(
QString(
"setEllipsoid: wrong format of parameter2 field: '%1'" ).arg( parameter2 ) );
220 QgsDebugMsg(
QString(
"setEllipsoid: a=%1, b=%2, 1/f=%3" ).arg( mSemiMajor ).arg( mSemiMinor ).arg( mInvFlattening ) );
224 QString proj4 =
"+proj=longlat +ellps=" + ellipsoid +
" +no_defs";
231 if ( destCRS.
srsid() == 0 )
234 .
arg(
QObject::tr(
"Generated CRS",
"A CRS automatically generated from layer info get this prefix for description" ),
256 mEllipsoid =
QString(
"PARAMETER:%1:%2" ).
arg( semiMajor ).
arg( semiMinor );
257 mSemiMajor = semiMajor;
258 mSemiMinor = semiMinor;
259 mInvFlattening = mSemiMajor / ( mSemiMajor - mSemiMinor );
274 if ( geomDimension <= 0 )
279 MeasureType measureType = type;
280 if ( measureType == Default )
282 measureType = ( geomDimension == 1 ? Length : Area );
285 if ( !mEllipsoidalMode || mEllipsoid ==
GEO_NONE )
288 if ( measureType == Length )
294 return geomV2->
area();
311 if ( measureType == Length )
362 return measure( geomV2, Area );
371 return measure( geomV2, Length );
380 if ( !geomV2 || geomV2->
dimension() < 2 )
385 if ( !mEllipsoidalMode || mEllipsoid ==
GEO_NONE )
403 surfaces.
append( static_cast<const QgsSurfaceV2*>( multiSurf->
geometryN( i ) ) );
409 for ( ; surfaceIt != surfaces.
constEnd(); ++surfaceIt )
420 length +=
measure( outerRing );
423 for (
int i = 0; i < nInnerRings; ++i )
441 curve->
points( linePointsV2 );
448 if ( points.
size() < 2 )
456 if ( mEllipsoidalMode && ( mEllipsoid !=
GEO_NONE ) )
457 p1 = mCoordTransform->
transform( points[0] );
463 if ( mEllipsoidalMode && ( mEllipsoid !=
GEO_NONE ) )
504 if ( mEllipsoidalMode && ( mEllipsoid !=
GEO_NONE ) )
563 for (
int idx = 0; idx < numRings; idx++ )
570 for (
int jdx = 0; jdx < nPoints; jdx++ )
576 wkbPtr +=
sizeof( double );
581 if ( mEllipsoidalMode && ( mEllipsoid !=
GEO_NONE ) )
588 if ( points.
size() > 2 )
634 curve->
points( linePointsV2 );
645 if ( mEllipsoidalMode && ( mEllipsoid !=
GEO_NONE ) )
673 if ( mEllipsoidalMode && ( mEllipsoid !=
GEO_NONE ) )
681 double dx = p2.
x() - p1.
x();
682 double dy = p2.
y() - p1.
y();
683 bearing = atan2( dx, dy );
695 double* course1,
double* course2 )
const
697 if ( p1.
x() == p2.
x() && p1.
y() == p2.
y() )
701 double a = mSemiMajor;
702 double b = mSemiMinor;
703 double f = 1 / mInvFlattening;
708 double L = p2_lon - p1_lon;
709 double U1 = atan(( 1 - f ) * tan( p1_lat ) );
710 double U2 = atan(( 1 - f ) * tan( p2_lat ) );
711 double sinU1 = sin( U1 ), cosU1 = cos( U1 );
712 double sinU2 = sin( U2 ), cosU2 = cos( U2 );
714 double lambdaP = 2 *
M_PI;
716 double sinLambda = 0;
717 double cosLambda = 0;
722 double cosSqAlpha = 0;
723 double cos2SigmaM = 0;
729 while ( qAbs( lambda - lambdaP ) > 1e-12 && --iterLimit > 0 )
731 sinLambda = sin( lambda );
732 cosLambda = cos( lambda );
733 tu1 = ( cosU2 * sinLambda );
734 tu2 = ( cosU1 * sinU2 - sinU1 * cosU2 * cosLambda );
735 sinSigma = sqrt( tu1 * tu1 + tu2 * tu2 );
736 cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
737 sigma = atan2( sinSigma, cosSigma );
738 alpha = asin( cosU1 * cosU2 * sinLambda / sinSigma );
739 cosSqAlpha = cos( alpha ) * cos( alpha );
740 cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
741 C = f / 16 * cosSqAlpha * ( 4 + f * ( 4 - 3 * cosSqAlpha ) );
743 lambda = L + ( 1 - C ) * f * sin( alpha ) *
744 ( sigma + C * sinSigma * ( cos2SigmaM + C * cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) ) );
747 if ( iterLimit == 0 )
750 double uSq = cosSqAlpha * ( a * a - b * b ) / ( b * b );
751 double A = 1 + uSq / 16384 * ( 4096 + uSq * ( -768 + uSq * ( 320 - 175 * uSq ) ) );
752 double B = uSq / 1024 * ( 256 + uSq * ( -128 + uSq * ( 74 - 47 * uSq ) ) );
753 double deltaSigma = B * sinSigma * ( cos2SigmaM + B / 4 * ( cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) -
754 B / 6 * cos2SigmaM * ( -3 + 4 * sinSigma * sinSigma ) * ( -3 + 4 * cos2SigmaM * cos2SigmaM ) ) );
755 double s = b * A * ( sigma - deltaSigma );
759 *course1 = atan2( tu1, tu2 );
764 *course2 = atan2( cosU1 * sinLambda, -sinU1 * cosU2 + cosU1 * sinU2 * cosLambda ) +
M_PI;
772 return sqrt(( p2.
x() - p1.
x() ) * ( p2.
x() - p1.
x() ) + ( p2.
y() - p1.
y() ) * ( p2.
y() - p1.
y() ) );
777 if ( points.
size() < 2 )
790 if ( mEllipsoidalMode && ( mEllipsoid !=
GEO_NONE ) )
819 double QgsDistanceArea::getQ(
double x )
const
826 return sinx *( 1 + sinx2 *( m_QA + sinx2 *( m_QB + sinx2 * m_QC ) ) );
830 double QgsDistanceArea::getQbar(
double x )
const
837 return cosx *( m_QbarA + cosx2 *( m_QbarB + cosx2 *( m_QbarC + cosx2 * m_QbarD ) ) );
849 double a2 = ( mSemiMajor * mSemiMajor );
850 double e2 = 1 - ( a2 / ( mSemiMinor * mSemiMinor ) );
858 m_AE = a2 * ( 1 - e2 );
860 m_QA = ( 2.0 / 3.0 ) * e2;
861 m_QB = ( 3.0 / 5.0 ) * e4;
862 m_QC = ( 4.0 / 7.0 ) * e6;
864 m_QbarA = -1.0 - ( 2.0 / 3.0 ) * e2 - ( 3.0 / 5.0 ) * e4 - ( 4.0 / 7.0 ) * e6;
865 m_QbarB = ( 2.0 / 9.0 ) * e2 + ( 2.0 / 5.0 ) * e4 + ( 4.0 / 7.0 ) * e6;
866 m_QbarC = - ( 3.0 / 25.0 ) * e4 - ( 12.0 / 35.0 ) * e6;
867 m_QbarD = ( 4.0 / 49.0 ) * e6;
869 m_Qp = getQ( M_PI / 2 );
870 m_E = 4 * M_PI * m_Qp * m_AE;
878 double x1, y1, x2, y2, dx, dy;
883 if (( ! mEllipsoidalMode ) || ( mEllipsoid ==
GEO_NONE ) )
887 int n = points.
size();
888 x2 =
DEG2RAD( points[n-1].x() );
889 y2 =
DEG2RAD( points[n-1].y() );
890 Qbar2 = getQbar( y2 );
894 for (
int i = 0; i < n; i++ )
902 Qbar2 = getQbar( y2 );
905 while ( x1 - x2 >
M_PI )
908 while ( x2 - x1 >
M_PI )
912 area += dx * ( m_Qp - getQ( y2 ) );
914 if (( dy = y2 - y1 ) != 0.0 )
915 area += dx * getQ( y2 ) - ( dx / dy ) * ( Qbar2 - Qbar1 );
917 if (( area *= m_AE ) < 0.0 )
927 if ( area > m_E / 2 )
939 size = points.
size();
942 for ( i = 0; i < size; i++ )
947 area = area + points[i].x() * points[( i+1 ) % size].y() - points[( i+1 ) % size].x() * points[i].y();
967 else if ( qAbs( value ) > 1000000.0 )
970 value = value / 1000000.0;
972 else if ( qAbs( value ) > 10000.0 )
975 value = value / 10000.0;
984 if ( keepBaseUnit || qAbs( value ) == 0.0 )
988 else if ( qAbs( value ) > 1000.0 )
991 value = value / 1000;
993 else if ( qAbs( value ) < 0.01 )
996 value = value * 1000;
998 else if ( qAbs( value ) < 0.1 )
1001 value = value * 100;
1012 if ( keepBaseUnit || qAbs( value ) <= 0.5*43560.0 )
1017 else if ( qAbs( value ) <= 0.5*5280.0*5280.0 )
1027 value /= 5280.0 * 5280.0;
1032 if ( qAbs( value ) <= 528.0 || keepBaseUnit )
1034 if ( qAbs( value ) == 1.0 )
1067 if ( qAbs( value ) == 1.0 )
1077 QgsDebugMsg(
QString(
"Error: not picked up map units - actual value = %1" ).arg( u ) );
1094 QgsDebugMsg(
"We're measuring on an ellipsoid or using projections, the system is returning meters" );
1096 else if ( mEllipsoidalMode && mEllipsoid ==
GEO_NONE )
1100 QgsDebugMsg(
"We're measuing on planimetric distance/area on given CRS, measured value is in CRS units" );
1106 factorUnits *= factorUnits;
1109 measure *= factorUnits;
1111 measureUnits = displayUnits;
double bearing(const QgsPoint &p1, const QgsPoint &p2) const
compute bearing - in radians
virtual QgsPolygonV2 * surfaceToPolygon() const =0
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.
long srsid() const
Get the SrsId - if possible.
UnitType
Map units that qgis supports.
QgsAbstractGeometryV2 * geometry() const
Returns the underlying geometry store.
virtual double length() const
Returns the length of the geometry.
QStringList split(const QString &sep, SplitBehavior behavior, Qt::CaseSensitivity cs) const
virtual void points(QList< QgsPointV2 > &pt) const =0
Returns a list of points within the curve.
void setSourceCrs(long srsid)
sets source spatial reference system (by QGIS CRS)
void computeAreaInit()
precalculates some values (must be called always when changing ellipsoid)
QgsCurveV2 * exteriorRing() const
Abstract base class for all geometries.
A geometry is the spatial representation of a feature.
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
double x() const
Get the x value of the point.
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".
bool createFromOgcWmsCrs(QString theCrs)
Set up this CRS from the given OGC CRS.
Q_DECL_DEPRECATED double measure(const QgsGeometry *geometry) const
General measurement (line distance or polygon area)
QString number(int n, int base)
bool createFromSrsId(const long theSrsId)
Set up this srs by fetching the appropriate information from the sqlite backend.
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
Line string geometry type.
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 ...
static void convertPointList(const QList< QgsPoint > &input, QList< QgsPointV2 > &output)
Upgrades a point list from QgsPoint to QgsPointV2.
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.
static void logMessage(const QString &message, const QString &tag=QString::null, MessageLevel level=WARNING)
add a message to the instance (and create it if necessary)
double measureLine(const QList< QgsPoint > &points) const
measures line
Multi surface geometry collection.
QString toString() const
String representation of the point (x,y)
double measurePerimeter(const QgsGeometry *geometry) const
measures perimeter of polygon
bool saveAsUserCRS(const QString &name)
Save the proj4-string as a custom CRS.
A class to represent a point.
double measureArea(const QgsGeometry *geometry) const
Measures the area of a geometry.
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.
int numGeometries() const
Returns the number of geometries within the collection.
double measureLength(const QgsGeometry *geometry) const
Measures the length of a geometry.
virtual QgsLineStringV2 * curveToLine() const =0
Returns a new line string geometry corresponding to a segmentized approximation of the curve...
QString mid(int position, int n) const
const QString & ellipsoid() const
returns ellipsoid's acronym
QgsPolygonV2 * surfaceToPolygon() const override
void setSourceAuthId(const QString &authid)
sets source spatial reference system by authid
Class for storing a coordinate reference system (CRS)
virtual int dimension() const =0
Returns the inherent dimension of the geometry.
virtual double area() const
Returns the area of the geometry.
static double fromUnitToUnitFactor(QGis::UnitType fromUnit, QGis::UnitType toUnit)
Returns the conversion factor between the specified units.
const QgsAbstractGeometryV2 * geometryN(int n) const
Returns a const reference to a geometry from within the collection.
QString left(int n) const
double y() const
Get the y value of the point.
static QString srsDbFilePath()
Returns the path to the srs.db file.
Custom exception class for Coordinate Reference System related exceptions.
const_iterator constEnd() const
static QString toLiteral(QGis::UnitType unit)
Provides the canonical name of the type value.
int numInteriorRings() const
virtual double perimeter() const
Returns the perimeter of the geometry.
const_iterator constBegin() const
Abstract base class for curved geometry type.
QString arg(qlonglong a, int fieldWidth, int base, const QChar &fillChar) const
QgsDistanceArea()
Constructor.
double computeDistanceFlat(const QgsPoint &p1, const QgsPoint &p2) const
uses flat / planimetric / Euclidean distance
QgsCurveV2 * interiorRing(int i) const
void setEllipsoidalMode(bool flag)
sets whether coordinates must be projected to ellipsoid before measuring
bool createFromProj4(const QString &theProjString)
Set up this srs by passing it a proj4 style formatted string.
QString toProj4() const
Get the Proj Proj4 string representation of this srs.
QByteArray toUtf8() const
QGis::UnitType mapUnits() const
Get the units that the projection is in.