35 #define M_PI 3.14159265358979323846
38 #define DEG2RAD(x) ((x)*M_PI/180)
65 if (
this == & origDA )
114 QString radius, parameter2;
120 sqlite3_stmt *myPreparedStatement;
134 if ( ellipsoid.startsWith(
"PARAMETER" ) )
136 QStringList paramList = ellipsoid.split(
":" );
137 bool semiMajorOk, semiMinorOk;
138 double semiMajor = paramList[1].toDouble( & semiMajorOk );
139 double semiMinor = paramList[2].toDouble( & semiMinorOk );
140 if ( semiMajorOk && semiMinorOk )
162 QString mySql =
"select radius, parameter2 from tbl_ellipsoid where acronym='" + ellipsoid +
"'";
163 myResult = sqlite3_prepare( myDatabase, mySql.toUtf8(), mySql.toUtf8().length(), &myPreparedStatement, &myTail );
165 if ( myResult == SQLITE_OK )
167 if ( sqlite3_step( myPreparedStatement ) == SQLITE_ROW )
169 radius = QString((
char * )sqlite3_column_text( myPreparedStatement, 0 ) );
170 parameter2 = QString((
char * )sqlite3_column_text( myPreparedStatement, 1 ) );
174 sqlite3_finalize( myPreparedStatement );
175 sqlite3_close( myDatabase );
178 if ( radius.isEmpty() || parameter2.isEmpty() )
180 QgsDebugMsg( QString(
"setEllipsoid: no row in tbl_ellipsoid for acronym '%1'" ).arg( ellipsoid ) );
185 if ( radius.left( 2 ) ==
"a=" )
189 QgsDebugMsg( QString(
"setEllipsoid: wrong format of radius field: '%1'" ).arg( radius ) );
196 if ( parameter2.left( 2 ) ==
"b=" )
201 else if ( parameter2.left( 3 ) ==
"rf=" )
208 QgsDebugMsg( QString(
"setEllipsoid: wrong format of parameter2 field: '%1'" ).arg( parameter2 ) );
216 QString proj4 =
"+proj=longlat +ellps=" + ellipsoid +
" +no_defs";
223 if ( destCRS.
srsid() == 0 )
225 QString myName = QString(
" * %1 (%2)" )
226 .arg(
QObject::tr(
"Generated CRS",
"A CRS automatically generated from layer info get this prefix for description" ) )
247 mEllipsoid = QString(
"PARAMETER:%1:%2" ).arg( semiMajor ).arg( semiMinor );
262 const unsigned char* wkb = geometry->
asWkb();
271 double res, resTotal = 0;
275 bool hasZptr =
false;
283 QgsDebugMsg(
"returning " + QString::number( res ) );
290 for ( i = 0; i < count; i++ )
295 QgsDebugMsg(
"returning " + QString::number( resTotal ) );
302 QgsDebugMsg(
"returning " + QString::number( res ) );
309 for ( i = 0; i < count; i++ )
319 QgsDebugMsg(
"returning " + QString::number( resTotal ) );
323 QgsDebugMsg( QString(
"measure: unexpected geometry type: %1" ).arg( wkbType ) );
333 const unsigned char* wkb = geometry->
asWkb();
341 double res = 0.0, resTotal = 0.0;
345 bool hasZptr =
false;
359 QgsDebugMsg(
"returning " + QString::number( res ) );
366 for ( i = 0; i < count; i++ )
376 QgsDebugMsg(
"returning " + QString::number( resTotal ) );
380 QgsDebugMsg( QString(
"measure: unexpected geometry type: %1" ).arg( wkbType ) );
392 QList<QgsPoint> points;
395 QgsDebugMsg(
"This feature WKB has " + QString::number( nPoints ) +
" points" );
397 for (
int i = 0; i < nPoints; ++i )
403 wkbPtr +=
sizeof( double );
415 if ( points.size() < 2 )
428 for ( QList<QgsPoint>::const_iterator i = points.begin(); i != points.end(); ++i )
471 QgsDebugMsgLevel( QString(
"New points are %1 and %2, calculating..." ).arg( pp1.
toString( 4 ) ).arg( pp2.toString( 4 ) ), 4 );
512 QList<QgsPoint> points;
522 for (
int idx = 0; idx < numRings; idx++ )
529 for (
int jdx = 0; jdx < nPoints; jdx++ )
535 wkbPtr +=
sizeof( double );
544 points.append( pnt );
547 if ( points.size() > 2 )
593 for ( QList<QgsPoint>::const_iterator i = points.begin(); i != points.end(); ++i )
626 double dx = p2.
x() - p1.
x();
627 double dy = p2.
y() - p1.
y();
628 bearing = atan2( dx, dy );
640 double* course1,
double* course2 )
642 if ( p1.
x() == p2.
x() && p1.
y() == p2.
y() )
653 double L = p2_lon - p1_lon;
654 double U1 = atan(( 1 - f ) * tan( p1_lat ) );
655 double U2 = atan(( 1 - f ) * tan( p2_lat ) );
656 double sinU1 = sin( U1 ), cosU1 = cos( U1 );
657 double sinU2 = sin( U2 ), cosU2 = cos( U2 );
659 double lambdaP = 2 *
M_PI;
661 double sinLambda = 0;
662 double cosLambda = 0;
667 double cosSqAlpha = 0;
668 double cos2SigmaM = 0;
674 while ( qAbs( lambda - lambdaP ) > 1e-12 && --iterLimit > 0 )
676 sinLambda = sin( lambda );
677 cosLambda = cos( lambda );
678 tu1 = ( cosU2 * sinLambda );
679 tu2 = ( cosU1 * sinU2 - sinU1 * cosU2 * cosLambda );
680 sinSigma = sqrt( tu1 * tu1 + tu2 * tu2 );
681 cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
682 sigma = atan2( sinSigma, cosSigma );
683 alpha = asin( cosU1 * cosU2 * sinLambda / sinSigma );
684 cosSqAlpha = cos( alpha ) * cos( alpha );
685 cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
686 C = f / 16 * cosSqAlpha * ( 4 + f * ( 4 - 3 * cosSqAlpha ) );
688 lambda = L + ( 1 - C ) * f * sin( alpha ) *
689 ( sigma + C * sinSigma * ( cos2SigmaM + C * cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) ) );
692 if ( iterLimit == 0 )
695 double uSq = cosSqAlpha * ( a * a - b * b ) / ( b * b );
696 double A = 1 + uSq / 16384 * ( 4096 + uSq * ( -768 + uSq * ( 320 - 175 * uSq ) ) );
697 double B = uSq / 1024 * ( 256 + uSq * ( -128 + uSq * ( 74 - 47 * uSq ) ) );
698 double deltaSigma = B * sinSigma * ( cos2SigmaM + B / 4 * ( cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) -
699 B / 6 * cos2SigmaM * ( -3 + 4 * sinSigma * sinSigma ) * ( -3 + 4 * cos2SigmaM * cos2SigmaM ) ) );
700 double s = b * A * ( sigma - deltaSigma );
704 *course1 = atan2( tu1, tu2 );
709 *course2 = atan2( cosU1 * sinLambda, -sinU1 * cosU2 + cosU1 * sinU2 * cosLambda ) +
M_PI;
717 return sqrt(( p2.
x() - p1.
x() ) * ( p2.
x() - p1.
x() ) + ( p2.
y() - p1.
y() ) * ( p2.
y() - p1.
y() ) );
722 if ( points.size() < 2 )
732 for ( QList<QgsPoint>::const_iterator i = points.begin(); i != points.end(); ++i )
771 return sinx *( 1 + sinx2 *(
m_QA + sinx2 *(
m_QB + sinx2 *
m_QC ) ) );
797 m_AE = a2 * ( 1 - e2 );
799 m_QA = ( 2.0 / 3.0 ) * e2;
800 m_QB = ( 3.0 / 5.0 ) * e4;
801 m_QC = ( 4.0 / 7.0 ) * e6;
803 m_QbarA = -1.0 - ( 2.0 / 3.0 ) * e2 - ( 3.0 / 5.0 ) * e4 - ( 4.0 / 7.0 ) * e6;
804 m_QbarB = ( 2.0 / 9.0 ) * e2 + ( 2.0 / 5.0 ) * e4 + ( 4.0 / 7.0 ) * e6;
805 m_QbarC = - ( 3.0 / 25.0 ) * e4 - ( 12.0 / 35.0 ) * e6;
817 double x1, y1, x2, y2, dx, dy;
826 int n = points.size();
827 x2 =
DEG2RAD( points[n-1].x() );
828 y2 =
DEG2RAD( points[n-1].y() );
833 for (
int i = 0; i < n; i++ )
844 while ( x1 - x2 >
M_PI )
847 while ( x2 - x1 >
M_PI )
853 if (( dy = y2 - y1 ) != 0.0 )
854 area += dx *
getQ( y2 ) - ( dx / dy ) * ( Qbar2 - Qbar1 );
856 if (( area *=
m_AE ) < 0.0 )
866 if ( area >
m_E / 2 )
878 size = points.size();
881 for ( i = 0; i <
size; i++ )
886 area = area + points[i].x() * points[( i+1 ) % size].y() - points[( i+1 ) % size].x() * points[i].y();
904 unitLabel = QObject::trUtf8(
" m²" );
906 else if ( qAbs( value ) > 1000000.0 )
908 unitLabel = QObject::trUtf8(
" km²" );
909 value = value / 1000000.0;
911 else if ( qAbs( value ) > 10000.0 )
914 value = value / 10000.0;
918 unitLabel = QObject::trUtf8(
" m²" );
923 if ( keepBaseUnit || qAbs( value ) == 0.0 )
927 else if ( qAbs( value ) > 1000.0 )
930 value = value / 1000;
932 else if ( qAbs( value ) < 0.01 )
935 value = value * 1000;
937 else if ( qAbs( value ) < 0.1 )
951 if ( keepBaseUnit || qAbs( value ) <= 0.5*43560.0 )
956 else if ( qAbs( value ) <= 0.5*5280.0*5280.0 )
966 value /= 5280.0 * 5280.0;
971 if ( qAbs( value ) <= 528.0 || keepBaseUnit )
973 if ( qAbs( value ) == 1.0 )
1006 if ( qAbs( value ) == 1.0 )
1015 QgsDebugMsg( QString(
"Error: not picked up map units - actual value = %1" ).arg( u ) );
1019 return QLocale::system().toString( value,
'f', decimals ) + unitLabel;
1033 QgsDebugMsg(
"We're measuring on an ellipsoid or using projections, the system is returning meters" );
1039 QgsDebugMsg(
"We're measuing on planimetric distance/area on given CRS, measured value is in CRS units" );
1045 factorUnits *= factorUnits;
1048 measure *= factorUnits;
1050 measureUnits = displayUnits;
double computePolygonFlatArea(const QList< QgsPoint > &points)
double computeDistance(const QList< QgsPoint > &points)
calculate distance with given coordinates (does not do a transform anymore)
~QgsDistanceArea()
Destructor.
bool saveAsUserCRS(QString name)
Copied from QgsCustomProjectionDialog /// Please refactor into SQL handler !!! ///.
QgsCoordinateTransform * mCoordTransform
used for transforming coordinates from source CRS to ellipsoid's coordinates
void setSourceCrs(long srsid)
sets source spatial reference system (by QGIS CRS)
void computeAreaInit()
precalculates some values (must be called always when changing ellipsoid)
WkbType
Used for symbology operations.
bool setEllipsoid(const QString &ellipsoid)
sets ellipsoid by its acronym
static void logMessage(QString message, QString tag=QString::null, MessageLevel level=WARNING)
add a message to the instance (and create it if necessary)
void _copy(const QgsDistanceArea &origDA)
Copy helper.
bool createFromOgcWmsCrs(QString theCrs)
Set up this CRS from the given OGC CRS.
bool createFromSrsId(const long theSrsId)
double measurePerimeter(QgsGeometry *geometry)
measures perimeter of polygon
double measurePolygon(const QList< QgsPoint > &points)
measures polygon area
#define QgsDebugMsgLevel(str, level)
double measure(QgsGeometry *geometry)
general measurement (line distance or polygon area)
double computePolygonArea(const QList< QgsPoint > &points)
calculates area of polygon on ellipsoid algorithm has been taken from GRASS: gis/area_poly1.c
double bearing(const QgsPoint &p1, const QgsPoint &p2)
compute bearing - in radians
const long GEOCRS_ID
Magic number for a geographic coord sys in QGIS srs.db tbl_srs.srs_id.
QString toString() const
String representation of the point (x,y)
A class to represent a point 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.
double mSemiMajor
ellipsoid parameters
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.
double measureLine(const QList< QgsPoint > &points)
measures line
const CORE_EXPORT QString GEO_NONE
Constant that holds the string representation for "No ellips/No CRS".
static double fromUnitToUnitFactor(QGis::UnitType fromUnit, QGis::UnitType toUnit)
Returns the conversion factor between the specified units.
UnitType
Map units that qgis supports.
void convertMeasurement(double &measure, QGis::UnitType &measureUnits, QGis::UnitType displayUnits, bool isArea)
Helper for conversion between physical units.
Custom exception class for Coordinate Reference System related exceptions.
bool mEllipsoidalMode
indicates whether we will transform coordinates
static QString toLiteral(QGis::UnitType unit)
Provides the canonical name of the type value.
double computeDistanceFlat(const QgsPoint &p1, const QgsPoint &p2)
uses flat / planimetric / Euclidean distance
QString mEllipsoid
ellipsoid acronym (from table tbl_ellipsoids)
const unsigned char * asWkb() const
Returns the buffer containing this geometry in WKB format.
QgsDistanceArea()
Constructor.
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.
double computeDistanceBearing(const QgsPoint &p1, const QgsPoint &p2, double *course1=NULL, double *course2=NULL)
calculates distance from two points on ellipsoid based on inverse Vincenty's formulae ...
void setSourceAuthId(QString authid)
sets source spatial reference system by authid
QGis::UnitType mapUnits() const