35 #define M_PI 3.14159265358979323846
38 #define DEG2RAD(x) ((x)*M_PI/180)
85 QString radius, parameter2;
91 sqlite3_stmt *myPreparedStatement;
95 if ( ellipsoid ==
"NONE" )
105 QgsDebugMsg( QString(
"Can't open database: %1" ).arg( sqlite3_errmsg( myDatabase ) ) );
111 QString mySql =
"select radius, parameter2 from tbl_ellipsoid where acronym='" + ellipsoid +
"'";
112 myResult = sqlite3_prepare( myDatabase, mySql.toUtf8(), mySql.toUtf8().length(), &myPreparedStatement, &myTail );
114 if ( myResult == SQLITE_OK )
116 if ( sqlite3_step( myPreparedStatement ) == SQLITE_ROW )
118 radius = QString((
char * )sqlite3_column_text( myPreparedStatement, 0 ) );
119 parameter2 = QString((
char * )sqlite3_column_text( myPreparedStatement, 1 ) );
123 sqlite3_finalize( myPreparedStatement );
124 sqlite3_close( myDatabase );
127 if ( radius.isEmpty() || parameter2.isEmpty() )
129 QgsDebugMsg( QString(
"setEllipsoid: no row in tbl_ellipsoid for acronym '%1'" ).arg( ellipsoid ) );
134 if ( radius.left( 2 ) ==
"a=" )
138 QgsDebugMsg( QString(
"setEllipsoid: wrong format of radius field: '%1'" ).arg( radius ) );
145 if ( parameter2.left( 2 ) ==
"b=" )
150 else if ( parameter2.left( 3 ) ==
"rf=" )
157 QgsDebugMsg( QString(
"setEllipsoid: wrong format of parameter2 field: '%1'" ).arg( parameter2 ) );
165 QString proj4 =
"+proj=longlat +ellps=" + ellipsoid +
" +no_defs";
185 unsigned char* wkb = geometry->
asWkb();
190 unsigned int wkbType;
191 double res, resTotal = 0;
194 memcpy( &wkbType, ( wkb + 1 ),
sizeof( wkbType ) );
197 bool hasZptr =
false;
205 QgsDebugMsg(
"returning " + QString::number( res ) );
211 count = *((
int* )( wkb + 5 ) );
213 for ( i = 0; i < count; i++ )
218 QgsDebugMsg(
"returning " + QString::number( resTotal ) );
225 QgsDebugMsg(
"returning " + QString::number( res ) );
231 count = *((
int* )( wkb + 5 ) );
233 for ( i = 0; i < count; i++ )
238 QgsDebugMsg(
"returning " + QString::number( resTotal ) );
242 QgsDebugMsg( QString(
"measure: unexpected geometry type: %1" ).arg( wkbType ) );
252 unsigned char* wkb = geometry->
asWkb();
257 unsigned int wkbType;
258 double res, resTotal = 0;
261 memcpy( &wkbType, ( wkb + 1 ),
sizeof( wkbType ) );
264 bool hasZptr =
false;
278 QgsDebugMsg(
"returning " + QString::number( res ) );
284 count = *((
int* )( wkb + 5 ) );
286 for ( i = 0; i < count; i++ )
291 QgsDebugMsg(
"returning " + QString::number( resTotal ) );
295 QgsDebugMsg( QString(
"measure: unexpected geometry type: %1" ).arg( wkbType ) );
303 unsigned char *ptr = feature + 5;
304 unsigned int nPoints = *((
int* )ptr );
307 QList<QgsPoint> points;
310 QgsDebugMsg(
"This feature WKB has " + QString::number( nPoints ) +
" points" );
312 for (
unsigned int i = 0; i < nPoints; ++i )
314 x = *((
double * ) ptr );
315 ptr +=
sizeof( double );
316 y = *((
double * ) ptr );
317 ptr +=
sizeof( double );
321 ptr +=
sizeof( double );
333 if ( points.size() < 2 )
346 for ( QList<QgsPoint>::const_iterator i = points.begin(); i != points.end(); ++i )
367 QgsLogger::warning( QObject::tr(
"Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
386 return sqrt(( p2.
x() - p1.
x() )*( p2.
x() - p1.
x() ) + ( p2.
y() - p1.
y() )*( p2.
y() - p1.
y() ) );
392 QgsLogger::warning( QObject::tr(
"Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
401 unsigned int numRings = *((
int* )( feature + 1 +
sizeof(
int ) ) );
407 unsigned char* ptr = feature + 1 + 2 *
sizeof( int );
409 QList<QgsPoint> points;
419 for (
unsigned int idx = 0; idx < numRings; idx++ )
421 int nPoints = *((
int* )ptr );
426 for (
int jdx = 0; jdx < nPoints; jdx++ )
428 x = *((
double * ) ptr );
429 ptr +=
sizeof( double );
430 y = *((
double * ) ptr );
431 ptr +=
sizeof( double );
435 ptr +=
sizeof( double );
444 points.append( pnt );
447 if ( points.size() > 2 )
480 QgsLogger::warning( QObject::tr(
"Caught a coordinate system exception while trying to transform a point. Unable to calculate polygon area or perimeter." ) );
495 for ( QList<QgsPoint>::const_iterator i = points.begin(); i != points.end(); ++i )
509 QgsLogger::warning( QObject::tr(
"Caught a coordinate system exception while trying to transform a point. Unable to calculate polygon area." ) );
528 double dx = p2.
x() - p1.
x();
529 double dy = p2.
y() - p1.
y();
530 bearing = atan2( dx, dy );
542 double* course1,
double* course2 )
544 if ( p1.
x() == p2.
x() && p1.
y() == p2.
y() )
555 double L = p2_lon - p1_lon;
556 double U1 = atan(( 1 - f ) * tan( p1_lat ) );
557 double U2 = atan(( 1 - f ) * tan( p2_lat ) );
558 double sinU1 = sin( U1 ), cosU1 = cos( U1 );
559 double sinU2 = sin( U2 ), cosU2 = cos( U2 );
561 double lambdaP = 2 *
M_PI;
563 double sinLambda = 0;
564 double cosLambda = 0;
569 double cosSqAlpha = 0;
570 double cos2SigmaM = 0;
576 while ( qAbs( lambda - lambdaP ) > 1e-12 && --iterLimit > 0 )
578 sinLambda = sin( lambda );
579 cosLambda = cos( lambda );
580 tu1 = ( cosU2 * sinLambda );
581 tu2 = ( cosU1 * sinU2 - sinU1 * cosU2 * cosLambda );
582 sinSigma = sqrt( tu1 * tu1 + tu2 * tu2 );
583 cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
584 sigma = atan2( sinSigma, cosSigma );
585 alpha = asin( cosU1 * cosU2 * sinLambda / sinSigma );
586 cosSqAlpha = cos( alpha ) * cos( alpha );
587 cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
588 C = f / 16 * cosSqAlpha * ( 4 + f * ( 4 - 3 * cosSqAlpha ) );
590 lambda = L + ( 1 - C ) * f * sin( alpha ) *
591 ( sigma + C * sinSigma * ( cos2SigmaM + C * cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) ) );
594 if ( iterLimit == 0 )
597 double uSq = cosSqAlpha * ( a * a - b * b ) / ( b * b );
598 double A = 1 + uSq / 16384 * ( 4096 + uSq * ( -768 + uSq * ( 320 - 175 * uSq ) ) );
599 double B = uSq / 1024 * ( 256 + uSq * ( -128 + uSq * ( 74 - 47 * uSq ) ) );
600 double deltaSigma = B * sinSigma * ( cos2SigmaM + B / 4 * ( cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) -
601 B / 6 * cos2SigmaM * ( -3 + 4 * sinSigma * sinSigma ) * ( -3 + 4 * cos2SigmaM * cos2SigmaM ) ) );
602 double s = b * A * ( sigma - deltaSigma );
606 *course1 = atan2( tu1, tu2 );
611 *course2 = atan2( cosU1 * sinLambda, -sinU1 * cosU2 + cosU1 * sinU2 * cosLambda ) +
M_PI;
631 return sinx *( 1 + sinx2 *(
m_QA + sinx2 *(
m_QB + sinx2 *
m_QC ) ) );
657 m_AE = a2 * ( 1 - e2 );
659 m_QA = ( 2.0 / 3.0 ) * e2;
660 m_QB = ( 3.0 / 5.0 ) * e4;
661 m_QC = ( 4.0 / 7.0 ) * e6;
663 m_QbarA = -1.0 - ( 2.0 / 3.0 ) * e2 - ( 3.0 / 5.0 ) * e4 - ( 4.0 / 7.0 ) * e6;
664 m_QbarB = ( 2.0 / 9.0 ) * e2 + ( 2.0 / 5.0 ) * e4 + ( 4.0 / 7.0 ) * e6;
665 m_QbarC = - ( 3.0 / 25.0 ) * e4 - ( 12.0 / 35.0 ) * e6;
676 double x1, y1, x2, y2, dx, dy;
685 int n = points.size();
686 x2 =
DEG2RAD( points[n-1].x() );
687 y2 =
DEG2RAD( points[n-1].y() );
692 for (
int i = 0; i < n; i++ )
703 while ( x1 - x2 >
M_PI )
706 while ( x2 - x1 >
M_PI )
712 if (( dy = y2 - y1 ) != 0.0 )
713 area += dx *
getQ( y2 ) - ( dx / dy ) * ( Qbar2 - Qbar1 );
715 if (( area *=
m_AE ) < 0.0 )
723 if ( area >
m_E ) area =
m_E;
724 if ( area >
m_E / 2 ) area =
m_E - area;
735 size = points.size();
738 for ( i = 0; i < size; i++ )
743 area = area + points[i].x() * points[( i+1 ) % size].y() - points[( i+1 ) % size].x() * points[i].y();
761 unitLabel = QObject::tr(
" m2" );
763 else if ( qAbs( value ) > 1000000.0 )
765 unitLabel = QObject::tr(
" km2" );
766 value = value / 1000000.0;
768 else if ( qAbs( value ) > 10000.0 )
770 unitLabel = QObject::tr(
" ha" );
771 value = value / 10000.0;
775 unitLabel = QObject::tr(
" m2" );
780 if ( keepBaseUnit || qAbs( value ) == 0.0 )
782 unitLabel = QObject::tr(
" m" );
784 else if ( qAbs( value ) > 1000.0 )
786 unitLabel = QObject::tr(
" km" );
787 value = value / 1000;
789 else if ( qAbs( value ) < 0.01 )
791 unitLabel = QObject::tr(
" mm" );
792 value = value * 1000;
794 else if ( qAbs( value ) < 0.1 )
796 unitLabel = QObject::tr(
" cm" );
801 unitLabel = QObject::tr(
" m" );
808 if ( keepBaseUnit || qAbs( value ) <= ( 528.0*528.0 ) )
810 unitLabel = QObject::tr(
" sq ft" );
814 unitLabel = QObject::tr(
" sq mile" );
815 value = value / ( 5280.0 * 5280.0 );
820 if ( qAbs( value ) <= 528.0 || keepBaseUnit )
822 if ( qAbs( value ) == 1.0 )
824 unitLabel = QObject::tr(
" foot" );
828 unitLabel = QObject::tr(
" feet" );
833 unitLabel = QObject::tr(
" mile" );
834 value = value / 5280.0;
841 unitLabel = QObject::tr(
" sq.deg." );
845 if ( qAbs( value ) == 1.0 )
846 unitLabel = QObject::tr(
" degree" );
848 unitLabel = QObject::tr(
" degrees" );
852 unitLabel = QObject::tr(
" unknown" );
854 QgsDebugMsg( QString(
"Error: not picked up map units - actual value = %1" ).arg( u ) );
858 return QLocale::system().toString( value,
'f', decimals ) + unitLabel;