35 #define M_PI 3.14159265358979323846
38 #define DEG2RAD(x) ((x)*M_PI/180)
65 if (
this == & origDA )
115 QString radius, parameter2;
121 sqlite3_stmt *myPreparedStatement;
135 if ( ellipsoid.startsWith(
"PARAMETER" ) )
137 QStringList paramList = ellipsoid.split(
":" );
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 );
179 if ( radius.isEmpty() || parameter2.isEmpty() )
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=" )
202 else if ( parameter2.left( 3 ) ==
"rf=" )
209 QgsDebugMsg( QString(
"setEllipsoid: wrong format of parameter2 field: '%1'" ).arg( parameter2 ) );
217 QString proj4 =
"+proj=longlat +ellps=" + ellipsoid +
" +no_defs";
224 if ( destCRS.
srsid() == 0 )
226 QString myName = QString(
" * %1 (%2)" )
227 .arg(
QObject::tr(
"Generated CRS",
"A CRS automatically generated from layer info get this prefix for description" ) )
248 mEllipsoid = QString(
"PARAMETER:%1:%2" ).arg( semiMajor ).arg( semiMinor );
263 const unsigned char* wkb = geometry->
asWkb();
272 double res, resTotal = 0;
276 bool hasZptr =
false;
284 QgsDebugMsg(
"returning " + QString::number( res ) );
291 for ( i = 0; i < count; i++ )
296 QgsDebugMsg(
"returning " + QString::number( resTotal ) );
303 QgsDebugMsg(
"returning " + QString::number( res ) );
310 for ( i = 0; i < count; i++ )
320 QgsDebugMsg(
"returning " + QString::number( resTotal ) );
324 QgsDebugMsg( QString(
"measure: unexpected geometry type: %1" ).arg( wkbType ) );
334 const unsigned char* wkb = geometry->
asWkb();
342 double res = 0.0, resTotal = 0.0;
346 bool hasZptr =
false;
360 QgsDebugMsg(
"returning " + QString::number( res ) );
367 for ( i = 0; i < count; i++ )
377 QgsDebugMsg(
"returning " + QString::number( resTotal ) );
381 QgsDebugMsg( QString(
"measure: unexpected geometry type: %1" ).arg( wkbType ) );
393 QList<QgsPoint> points;
396 QgsDebugMsg(
"This feature WKB has " + QString::number( nPoints ) +
" points" );
398 for (
int i = 0; i < nPoints; ++i )
404 wkbPtr +=
sizeof( double );
416 if ( points.size() < 2 )
429 for ( QList<QgsPoint>::const_iterator i = points.begin(); i != points.end(); ++i )
472 QgsDebugMsgLevel( QString(
"New points are %1 and %2, calculating..." ).arg( pp1.
toString( 4 ) ).arg( pp2.toString( 4 ) ), 4 );
478 result = sqrt(( p2.
x() - p1.
x() ) * ( p2.
x() - p1.
x() ) + ( p2.
y() - p1.
y() ) * ( p2.
y() - p1.
y() ) );
513 QList<QgsPoint> points;
523 for (
int idx = 0; idx < numRings; idx++ )
530 for (
int jdx = 0; jdx < nPoints; jdx++ )
536 wkbPtr +=
sizeof( double );
545 points.append( pnt );
548 if ( points.size() > 2 )
594 for ( QList<QgsPoint>::const_iterator i = points.begin(); i != points.end(); ++i )
627 double dx = p2.
x() - p1.
x();
628 double dy = p2.
y() - p1.
y();
629 bearing = atan2( dx, dy );
641 double* course1,
double* course2 )
643 if ( p1.
x() == p2.
x() && p1.
y() == p2.
y() )
654 double L = p2_lon - p1_lon;
655 double U1 = atan(( 1 - f ) * tan( p1_lat ) );
656 double U2 = atan(( 1 - f ) * tan( p2_lat ) );
657 double sinU1 = sin( U1 ), cosU1 = cos( U1 );
658 double sinU2 = sin( U2 ), cosU2 = cos( U2 );
660 double lambdaP = 2 *
M_PI;
662 double sinLambda = 0;
663 double cosLambda = 0;
668 double cosSqAlpha = 0;
669 double cos2SigmaM = 0;
675 while ( qAbs( lambda - lambdaP ) > 1e-12 && --iterLimit > 0 )
677 sinLambda = sin( lambda );
678 cosLambda = cos( lambda );
679 tu1 = ( cosU2 * sinLambda );
680 tu2 = ( cosU1 * sinU2 - sinU1 * cosU2 * cosLambda );
681 sinSigma = sqrt( tu1 * tu1 + tu2 * tu2 );
682 cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
683 sigma = atan2( sinSigma, cosSigma );
684 alpha = asin( cosU1 * cosU2 * sinLambda / sinSigma );
685 cosSqAlpha = cos( alpha ) * cos( alpha );
686 cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
687 C = f / 16 * cosSqAlpha * ( 4 + f * ( 4 - 3 * cosSqAlpha ) );
689 lambda = L + ( 1 - C ) * f * sin( alpha ) *
690 ( sigma + C * sinSigma * ( cos2SigmaM + C * cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) ) );
693 if ( iterLimit == 0 )
696 double uSq = cosSqAlpha * ( a * a - b * b ) / ( b * b );
697 double A = 1 + uSq / 16384 * ( 4096 + uSq * ( -768 + uSq * ( 320 - 175 * uSq ) ) );
698 double B = uSq / 1024 * ( 256 + uSq * ( -128 + uSq * ( 74 - 47 * uSq ) ) );
699 double deltaSigma = B * sinSigma * ( cos2SigmaM + B / 4 * ( cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) -
700 B / 6 * cos2SigmaM * ( -3 + 4 * sinSigma * sinSigma ) * ( -3 + 4 * cos2SigmaM * cos2SigmaM ) ) );
701 double s = b * A * ( sigma - deltaSigma );
705 *course1 = atan2( tu1, tu2 );
710 *course2 = atan2( cosU1 * sinLambda, -sinU1 * cosU2 + cosU1 * sinU2 * cosLambda ) +
M_PI;
730 return sinx *( 1 + sinx2 *(
m_QA + sinx2 *(
m_QB + sinx2 *
m_QC ) ) );
756 m_AE = a2 * ( 1 - e2 );
758 m_QA = ( 2.0 / 3.0 ) * e2;
759 m_QB = ( 3.0 / 5.0 ) * e4;
760 m_QC = ( 4.0 / 7.0 ) * e6;
762 m_QbarA = -1.0 - ( 2.0 / 3.0 ) * e2 - ( 3.0 / 5.0 ) * e4 - ( 4.0 / 7.0 ) * e6;
763 m_QbarB = ( 2.0 / 9.0 ) * e2 + ( 2.0 / 5.0 ) * e4 + ( 4.0 / 7.0 ) * e6;
764 m_QbarC = - ( 3.0 / 25.0 ) * e4 - ( 12.0 / 35.0 ) * e6;
776 double x1, y1, x2, y2, dx, dy;
785 int n = points.size();
786 x2 =
DEG2RAD( points[n-1].x() );
787 y2 =
DEG2RAD( points[n-1].y() );
792 for (
int i = 0; i < n; i++ )
803 while ( x1 - x2 >
M_PI )
806 while ( x2 - x1 >
M_PI )
812 if (( dy = y2 - y1 ) != 0.0 )
813 area += dx *
getQ( y2 ) - ( dx / dy ) * ( Qbar2 - Qbar1 );
815 if (( area *=
m_AE ) < 0.0 )
825 if ( area >
m_E / 2 )
837 size = points.size();
840 for ( i = 0; i <
size; i++ )
845 area = area + points[i].x() * points[( i+1 ) % size].y() - points[( i+1 ) % size].x() * points[i].y();
863 unitLabel = QObject::trUtf8(
" m²" );
865 else if ( qAbs( value ) > 1000000.0 )
867 unitLabel = QObject::trUtf8(
" km²" );
868 value = value / 1000000.0;
870 else if ( qAbs( value ) > 10000.0 )
873 value = value / 10000.0;
877 unitLabel = QObject::trUtf8(
" m²" );
882 if ( keepBaseUnit || qAbs( value ) == 0.0 )
886 else if ( qAbs( value ) > 1000.0 )
889 value = value / 1000;
891 else if ( qAbs( value ) < 0.01 )
894 value = value * 1000;
896 else if ( qAbs( value ) < 0.1 )
910 if ( keepBaseUnit || qAbs( value ) <= 0.5*43560.0 )
915 else if ( qAbs( value ) <= 0.5*5280.0*5280.0 )
925 value /= 5280.0 * 5280.0;
930 if ( qAbs( value ) <= 528.0 || keepBaseUnit )
932 if ( qAbs( value ) == 1.0 )
965 if ( qAbs( value ) == 1.0 )
974 QgsDebugMsg( QString(
"Error: not picked up map units - actual value = %1" ).arg( u ) );
978 return QLocale::system().toString( value,
'f', decimals ) + unitLabel;
992 QgsDebugMsg(
"We're measuring on an ellipsoid or using projections, the system is returning meters" );
998 QgsDebugMsg(
"We're measuing on planimetric distance/area on given CRS, measured value is in CRS units" );
1004 factorUnits *= factorUnits;
1007 measure *= factorUnits;
1009 measureUnits = displayUnits;
double computePolygonFlatArea(const QList< QgsPoint > &points)
~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.
long mSourceRefSys
id of the source spatial reference system
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