1 /***************************************************************************
2 qgswcsprovider.cpp - QGIS Data provider for
3 OGC Web Coverage Service layers
4 -------------------
5 begin : 2 July, 2012
6 copyright : (C) 2012 by Radim Blazek
7 email : radim dot blazek at gmail.com
8
9 Based on qgswmsprovider.cpp written by Brendan Morley.
10
11 ***************************************************************************/
12
13 /***************************************************************************
14 * *
15 * This program is free software; you can redistribute it and/or modify *
16 * it under the terms of the GNU General Public License as published by *
17 * the Free Software Foundation; either version 2 of the License, or *
18 * (at your option) any later version. *
19 * *
20 ***************************************************************************/
21
22 #include "qgslogger.h"
23 #include "qgswcsprovider.h"
24 #include "qgscoordinatetransform.h"
25 #include "qgsdatasourceuri.h"
26 #include "qgsrasteridentifyresult.h"
27 #include "qgsrasterlayer.h"
28 #include "qgsrectangle.h"
29 #include "qgscoordinatereferencesystem.h"
30 #include "qgsnetworkaccessmanager.h"
31 #include "qgsnetworkreplyparser.h"
32 #include "qgsmessagelog.h"
33 #include "qgsexception.h"
34 #include "qgswcsdataitems.h"
35
36 #include <QNetworkRequest>
37 #include <QNetworkReply>
38 #include <QNetworkProxy>
39
40 #include <QUrl>
41 #include <QEventLoop>
42 #include <QFile>
43
44 #ifdef QGISDEBUG
45 #include <QDir>
46 #endif
47
48 #include "gdalwarper.h"
49 #include "ogr_srs_api.h"
50 #include "cpl_conv.h"
51 #include "cpl_string.h"
52
53 #define ERR(message) QGS_ERROR_MESSAGE(message,"WCS provider")
54 #define SRVERR(message) QGS_ERROR_MESSAGE(message,"WCS server")
55 #define QGS_ERROR(message) QgsError(message,"WCS provider")
56
57 QString QgsWcsProvider::WCS_KEY = QStringLiteral( "wcs" );
58 QString QgsWcsProvider::WCS_DESCRIPTION = QStringLiteral( "OGC Web Coverage Service version 1.0/1.1 data provider" );
59
60 static QString DEFAULT_LATLON_CRS = QStringLiteral( "CRS:84" );
61
62 // TODO: colortable - use common baseclass with gdal, mapserver does not support http://trac.osgeo.org/mapserver/ticket/1671
63
QgsWcsProvider(const QString & uri,const ProviderOptions & options,QgsDataProvider::ReadFlags flags)64 QgsWcsProvider::QgsWcsProvider( const QString &uri, const ProviderOptions &options, QgsDataProvider::ReadFlags flags )
65 : QgsRasterDataProvider( uri, options, flags )
66 , mCachedViewExtent( 0 )
67 {
68 QgsDebugMsg( "constructing with uri '" + mHttpUri + "'." );
69
70 mValid = false;
71 mCachedMemFilename = QStringLiteral( "/vsimem/qgis/wcs/%0.dat" ).arg( reinterpret_cast<std::uintptr_t>( this ) );
72
73 if ( !parseUri( uri ) ) return;
74
75 // GetCapabilities and DescribeCoverage
76 // TODO(?): do only DescribeCoverage to avoid one request
77 // We need to get at least server version, which is not set in of URI (if not part of url)
78 // and probably also rangeSet
79
80 QgsDataSourceUri capabilitiesUri;
81 capabilitiesUri.setEncodedUri( uri );
82 // remove non relevant params
83 capabilitiesUri.removeParam( QStringLiteral( "identifier" ) );
84 capabilitiesUri.removeParam( QStringLiteral( "crs" ) );
85 capabilitiesUri.removeParam( QStringLiteral( "format" ) );
86 // TODO: check if successful (add return to capabilities)
87 mCapabilities.setUri( capabilitiesUri );
88
89 // 1.0 get additional coverage info
90 if ( !mCapabilities.describeCoverage( mIdentifier ) )
91 {
92 appendError( ERR( tr( "Cannot describe coverage" ) ) );
93 return;
94 }
95
96 mCoverageSummary = mCapabilities.coverage( mIdentifier );
97 if ( !mCoverageSummary.valid )
98 {
99 // Should not happen if describeCoverage() did not fail
100 appendError( ERR( tr( "Coverage not found" ) ) );
101 return;
102 }
103
104 // It may happen that format is empty (e.g. uri created in python script),
105 // in that casei select one from available formats
106 if ( mFormat.isEmpty() )
107 {
108 // TIFF is known by GDAL
109 mFormat = mCoverageSummary.supportedFormat.filter( QStringLiteral( "tif" ), Qt::CaseInsensitive ).value( 0 );
110 }
111 if ( mFormat.isEmpty() )
112 {
113 // Take the first if TIFF was not found
114 mFormat = mCoverageSummary.supportedFormat.value( 0 );
115 }
116
117 // We cannot continue without format, it is required
118 if ( mFormat.isEmpty() ) return;
119
120 // It could happen (usually not with current QgsWCSSourceSelect if at least
121 // one CRS is available) that crs is not set in uri, in that case we
122 // use the native, if available or WGS84 or the first supported
123 if ( mCoverageCrs.isEmpty() )
124 {
125 QgsDebugMsg( "nativeCrs = " + mCoverageSummary.nativeCrs );
126 QgsDebugMsg( "supportedCrs = " + mCoverageSummary.supportedCrs.join( "," ) );
127 if ( !mCoverageSummary.nativeCrs.isEmpty() )
128 {
129 setCoverageCrs( mCoverageSummary.nativeCrs );
130 }
131 else if ( mCoverageSummary.supportedCrs.contains( QStringLiteral( "EPSG:4326" ), Qt::CaseInsensitive ) )
132 {
133 setCoverageCrs( QStringLiteral( "EPSG:4326" ) );
134 }
135 else if ( !mCoverageSummary.supportedCrs.isEmpty() )
136 {
137 setCoverageCrs( mCoverageSummary.supportedCrs.value( 0 ) );
138 }
139 }
140 QgsDebugMsg( "mCoverageCrs = " + mCoverageCrs );
141
142 // It may happen that coverage CRS is not given or it is unknown
143 // in that case we continue without CRS and user is asked for it
144 //if ( mCoverageCrs.isEmpty() ) return;
145
146 // Native size
147 mWidth = mCoverageSummary.width;
148 mHeight = mCoverageSummary.height;
149 mHasSize = mCoverageSummary.hasSize;
150
151 QgsDebugMsg( QStringLiteral( "mWidth = %1 mHeight = %2" ).arg( mWidth ).arg( mHeight ) );
152
153 // TODO: Consider if/how to recalculate mWidth, mHeight if non native CRS is used
154
155 if ( !calculateExtent() )
156 {
157 appendError( ERR( tr( "Cannot calculate extent" ) ) );
158 return;
159 }
160
161 // Get small piece of coverage to find GDAL data type and number of bands
162 int bandNo = 0; // All bands
163 int width;
164 int height;
165 QString crs;
166 QgsRectangle box; // box to use to calc extent
167 // Prefer native CRS
168 if ( !mCoverageSummary.nativeCrs.isEmpty() &&
169 !mCoverageSummary.nativeBoundingBox.isEmpty() &&
170 mCoverageSummary.supportedCrs.contains( mCoverageSummary.nativeCrs ) &&
171 mHasSize )
172 {
173 box = mCoverageSummary.nativeBoundingBox;
174 width = mWidth;
175 height = mHeight;
176 crs = mCoverageSummary.nativeCrs;
177 }
178 else
179 {
180 box = mCoverageExtent;
181 if ( mHasSize )
182 {
183 width = mWidth;
184 height = mHeight;
185 }
186 else
187 {
188 // just a number to get smaller piece of coverage
189 width = 1000;
190 height = 1000;
191 }
192 }
193 double xRes = box.width() / width;
194 double yRes = box.height() / height;
195 QgsPointXY p = box.center();
196
197 // width and height different to recognize rotation
198 int requestWidth = 6;
199 int requestHeight = 3;
200
201 // extent to be used for test request
202 double halfWidth = xRes * ( requestWidth / 2. );
203 double halfHeight = yRes * ( requestHeight / 2. );
204
205 // Using non native CRS (if we don't know which is native) it could easily happen,
206 // that a small part of bbox in request CRS near margin falls outside
207 // coverage native bbox and server reports error => take a piece from center
208 QgsRectangle extent = QgsRectangle( p.x() - halfWidth, p.y() - halfHeight, p.x() + halfWidth, p.y() + halfHeight );
209
210 getCache( bandNo, extent, requestWidth, requestHeight, crs );
211
212 if ( !mCachedGdalDataset )
213 {
214 setError( mCachedError );
215 appendError( ERR( tr( "Cannot get test dataset." ) ) );
216 return;
217 }
218
219 mBandCount = GDALGetRasterCount( mCachedGdalDataset.get() );
220 QgsDebugMsg( QStringLiteral( "mBandCount = %1" ).arg( mBandCount ) );
221
222 // Check for server particularities (bbox, rotation)
223 int responseWidth = GDALGetRasterXSize( mCachedGdalDataset.get() );
224 int responseHeight = GDALGetRasterYSize( mCachedGdalDataset.get() );
225
226 QgsDebugMsg( QStringLiteral( "requestWidth = %1 requestHeight = %2 responseWidth = %3 responseHeight = %4)" ).arg( requestWidth ).arg( requestHeight ).arg( responseWidth ).arg( responseHeight ) );
227 // GeoServer and ArcGIS are using for 1.1 box "pixel" edges
228 // Mapserver is using pixel centers according to 1.1. specification
229 if ( ( responseWidth == requestWidth - 1 && responseHeight == requestHeight - 1 ) ||
230 ( responseWidth == requestHeight - 1 && responseHeight == requestWidth - 1 ) )
231 {
232 mFixBox = true;
233 QgsDebugMsg( QStringLiteral( "Test response size is smaller by pixel, using mFixBox" ) );
234 }
235 // Geoserver is giving rotated raster for geographic CRS - switched axis,
236 // Geoserver developers argue that changed axis order applies also to
237 // returned raster, that is exaggerated IMO but we have to handle that.
238 if ( ( responseWidth == requestHeight && responseHeight == requestWidth ) ||
239 ( responseWidth == requestHeight - 1 && responseHeight == requestWidth - 1 ) )
240 {
241 mFixRotate = true;
242 QgsDebugMsg( QStringLiteral( "Test response is rotated, using mFixRotate" ) );
243 }
244
245 // Get types
246 // TODO: we are using the same data types like GDAL (not wider like GDAL provider)
247 // with expectation to replace 'no data' values by NaN
248 mSrcGdalDataType.reserve( mBandCount );
249 mGdalDataType.reserve( mBandCount );
250 for ( int i = 1; i <= mBandCount; i++ )
251 {
252 GDALRasterBandH gdalBand = GDALGetRasterBand( mCachedGdalDataset.get(), i );
253 GDALDataType myGdalDataType = GDALGetRasterDataType( gdalBand );
254
255 QgsDebugMsg( QStringLiteral( "myGdalDataType[%1] = %2" ).arg( i - 1 ).arg( myGdalDataType ) );
256 mSrcGdalDataType.append( myGdalDataType );
257
258 // UMN Mapserver does not automatically set null value, METADATA wcs_rangeset_nullvalue must be used
259 // http://lists.osgeo.org/pipermail/mapserver-users/2010-April/065328.html
260
261 // TODO: This could be shared with GDAL provider
262 int isValid = false;
263 double myNoDataValue = GDALGetRasterNoDataValue( gdalBand, &isValid );
264 if ( isValid )
265 {
266 QgsDebugMsg( QStringLiteral( "GDALGetRasterNoDataValue = %1" ).arg( myNoDataValue ) );
267 myNoDataValue = QgsRaster::representableValue( myNoDataValue, dataTypeFromGdal( myGdalDataType ) );
268 mSrcNoDataValue.append( myNoDataValue );
269 mSrcHasNoDataValue.append( true );
270 mUseSrcNoDataValue.append( true );
271 }
272 else
273 {
274 mSrcNoDataValue.append( std::numeric_limits<double>::quiet_NaN() );
275 mSrcHasNoDataValue.append( false );
276 mUseSrcNoDataValue.append( false );
277 }
278 // It may happen that nodata value given by GDAL is wrong and it has to be
279 // disabled by user, in that case we need another value to be used for nodata
280 // (for reprojection for example) -> always internally represent as wider type
281 // with mInternalNoDataValue in reserve.
282 // No retyping, no internal values for now
283 #if 0
284 int myInternalGdalDataType = myGdalDataType;
285 double myInternalNoDataValue;
286 switch ( srcDataType( i ) )
287 {
288 case Qgis::Byte:
289 myInternalNoDataValue = -32768.0;
290 myInternalGdalDataType = GDT_Int16;
291 break;
292 case Qgis::Int16:
293 myInternalNoDataValue = -2147483648.0;
294 myInternalGdalDataType = GDT_Int32;
295 break;
296 case Qgis::UInt16:
297 myInternalNoDataValue = -2147483648.0;
298 myInternalGdalDataType = GDT_Int32;
299 break;
300 case Qgis::Int32:
301 // We believe that such values is no used in real data
302 myInternalNoDataValue = -2147483648.0;
303 break;
304 case Qgis::UInt32:
305 // We believe that such values is no used in real data
306 myInternalNoDataValue = 4294967295.0;
307 break;
308 default: // Float32, Float64
309 //myNoDataValue = std::numeric_limits<int>::max();
310 // NaN should work well
311 myInternalNoDataValue = std::numeric_limits<double>::quiet_NaN();
312 }
313 mGdalDataType.append( myInternalGdalDataType );
314 mInternalNoDataValue.append( myInternalNoDataValue );
315 QgsDebugMsg( QStringLiteral( "mInternalNoDataValue[%1] = %2" ).arg( i - 1 ).arg( mInternalNoDataValue[i - 1] ) );
316 #endif
317 mGdalDataType.append( myGdalDataType );
318
319 #if 0
320 // TODO: what to do if null values from DescribeCoverage differ?
321 if ( !mCoverageSummary.nullValues.contains( myNoDataValue ) )
322 {
323 QgsDebugMsg( QStringLiteral( "noDataValue %1 is missing in nullValues from CoverageDescription" ).arg( myNoDataValue ) );
324 }
325 #endif
326
327 QgsDebugMsg( QStringLiteral( "mSrcGdalDataType[%1] = %2" ).arg( i - 1 ).arg( mSrcGdalDataType.at( i - 1 ) ) );
328 QgsDebugMsg( QStringLiteral( "mGdalDataType[%1] = %2" ).arg( i - 1 ).arg( mGdalDataType.at( i - 1 ) ) );
329 QgsDebugMsg( QStringLiteral( "mSrcNoDataValue[%1] = %2" ).arg( i - 1 ).arg( mSrcNoDataValue.at( i - 1 ) ) );
330
331 // Create and store color table
332 // TODO: never tested because mapserver (6.0.3) does not support color tables
333 mColorTables.append( QgsGdalProviderBase::colorTable( mCachedGdalDataset.get(), i ) );
334 }
335
336 clearCache();
337
338 // Block size is used for for statistics
339 // TODO: How to find maximum block size supported by server?
340 if ( mHasSize )
341 {
342 // This is taken from GDAL, how they come to these numbers?
343 if ( mWidth > 1800 ) mXBlockSize = 1024;
344 else mXBlockSize = mWidth;
345
346 if ( mHeight > 900 ) mYBlockSize = 512;
347 else mYBlockSize = mHeight;
348 }
349
350 mValid = true;
351 QgsDebugMsg( QStringLiteral( "Constructed OK, provider valid." ) );
352 }
353
QgsWcsProvider(const QgsWcsProvider & other,const QgsDataProvider::ProviderOptions & providerOptions)354 QgsWcsProvider::QgsWcsProvider( const QgsWcsProvider &other, const QgsDataProvider::ProviderOptions &providerOptions )
355 : QgsRasterDataProvider( other.dataSourceUri(), providerOptions )
356 , mHttpUri( other.mHttpUri )
357 , mBaseUrl( other.mBaseUrl )
358 , mIdentifier( other.mIdentifier )
359 , mTime( other.mTime )
360 , mFormat( other.mFormat )
361 , mValid( other.mValid )
362 , mCapabilities( other.mCapabilities )
363 , mCoverageSummary( other.mCoverageSummary )
364 , mSrid( other.mSrid )
365 , mCoverageExtent( other.mCoverageExtent )
366 , mWidth( other.mWidth )
367 , mHeight( other.mHeight )
368 , mXBlockSize( other.mXBlockSize )
369 , mYBlockSize( other.mYBlockSize )
370 , mHasSize( other.mHasSize )
371 , mBandCount( other.mBandCount )
372 , mGdalDataType( other.mGdalDataType )
373 , mSrcGdalDataType( other.mSrcGdalDataType )
374 , mColorTables( other.mColorTables )
375 , mExtentForLayer( other.mExtentForLayer )
376 , mCrsForLayer( other.mCrsForLayer )
377 , mQueryableForLayer( other.mQueryableForLayer )
378 , mCoverageCrs( other.mCoverageCrs )
379 // intentionally omitted:
380 // - mCachedData
381 // - mCachedMemFilename
382 // - mCachedMemFile
383 // - mCachedGdalDataset
384 // - mCachedError
385 // - mCachedViewExtent
386 // - mCachedViewWidth
387 // - mCachedViewHeight
388 , mMaxWidth( other.mMaxWidth )
389 , mMaxHeight( other.mMaxHeight )
390 // intentionally omitted:
391 // - mErrorCaption
392 // - mError
393 // - mErrorFormat
394 , mCoordinateTransform( other.mCoordinateTransform )
395 , mExtentDirty( other.mExtentDirty )
396 , mGetFeatureInfoUrlBase( other.mGetFeatureInfoUrlBase )
397 , mServiceMetadataURL( other.mServiceMetadataURL )
398 , mAuth( other.mAuth )
399 , mIgnoreGetCoverageUrl( other.mIgnoreGetCoverageUrl )
400 , mIgnoreAxisOrientation( other.mIgnoreAxisOrientation )
401 , mInvertAxisOrientation( other.mInvertAxisOrientation )
402 , mCrs( other.mCrs )
403 , mFixBox( other.mFixBox )
404 , mFixRotate( other.mFixRotate )
405 , mCacheLoadControl( other.mCacheLoadControl )
406 {
407 mCachedMemFilename = QStringLiteral( "/vsimem/qgis/wcs/%0.dat" ).arg( reinterpret_cast<std::uintptr_t>( this ) );
408 }
409
410
parseUri(const QString & uriString)411 bool QgsWcsProvider::parseUri( const QString &uriString )
412 {
413
414 QgsDebugMsg( "uriString = " + uriString );
415 QgsDataSourceUri uri;
416 uri.setEncodedUri( uriString );
417
418 mMaxWidth = 0;
419 mMaxHeight = 0;
420
421 mHttpUri = uri.param( QStringLiteral( "url" ) );
422 mBaseUrl = prepareUri( mHttpUri );
423 QgsDebugMsg( "mBaseUrl = " + mBaseUrl );
424
425 mIgnoreGetCoverageUrl = uri.hasParam( QStringLiteral( "IgnoreGetMapUrl" ) );
426 mIgnoreAxisOrientation = uri.hasParam( QStringLiteral( "IgnoreAxisOrientation" ) ); // must be before parsing!
427 mInvertAxisOrientation = uri.hasParam( QStringLiteral( "InvertAxisOrientation" ) ); // must be before parsing!
428
429 mAuth.mUserName = uri.username();
430 QgsDebugMsg( "set username to " + mAuth.mUserName );
431
432 mAuth.mPassword = uri.password();
433 QgsDebugMsg( "set password to " + mAuth.mPassword );
434
435 if ( !uri.authConfigId().isEmpty() )
436 {
437 mAuth.mAuthCfg = uri.authConfigId();
438 }
439 QgsDebugMsg( "set authcfg to " + mAuth.mAuthCfg );
440
441 mIdentifier = uri.param( QStringLiteral( "identifier" ) );
442
443 mTime = uri.param( QStringLiteral( "time" ) );
444
445 setFormat( uri.param( QStringLiteral( "format" ) ) );
446
447 if ( !uri.param( QStringLiteral( "crs" ) ).isEmpty() )
448 {
449 setCoverageCrs( uri.param( QStringLiteral( "crs" ) ) );
450 }
451
452 QString cache = uri.param( QStringLiteral( "cache" ) );
453 if ( !cache.isEmpty() )
454 {
455 mCacheLoadControl = QgsNetworkAccessManager::cacheLoadControlFromName( cache );
456 }
457 QgsDebugMsg( QStringLiteral( "mCacheLoadControl = %1" ).arg( mCacheLoadControl ) );
458
459 return true;
460 }
461
prepareUri(QString uri) const462 QString QgsWcsProvider::prepareUri( QString uri ) const
463 {
464 if ( !uri.contains( '?' ) )
465 {
466 uri.append( '?' );
467 }
468 else if ( uri.right( 1 ) != QLatin1String( "?" ) && uri.right( 1 ) != QLatin1String( "&" ) )
469 {
470 uri.append( '&' );
471 }
472
473 return uri;
474 }
475
~QgsWcsProvider()476 QgsWcsProvider::~QgsWcsProvider()
477 {
478 QgsDebugMsg( QStringLiteral( "deconstructing." ) );
479
480 // Dispose of any cached image as created by draw()
481 clearCache();
482 }
483
clone() const484 QgsWcsProvider *QgsWcsProvider::clone() const
485 {
486 QgsDataProvider::ProviderOptions providerOptions;
487 providerOptions.transformContext = transformContext();
488 QgsWcsProvider *provider = new QgsWcsProvider( *this, providerOptions );
489 provider->copyBaseSettings( *this );
490 return provider;
491 }
492
baseUrl() const493 QString QgsWcsProvider::baseUrl() const
494 {
495 return mBaseUrl;
496 }
497
format() const498 QString QgsWcsProvider::format() const
499 {
500 return mFormat;
501 }
502
setFormat(QString const & format)503 void QgsWcsProvider::setFormat( QString const &format )
504 {
505 QgsDebugMsg( "Setting format to " + format + '.' );
506 mFormat = format;
507 }
508
509
setCoverageCrs(QString const & crs)510 void QgsWcsProvider::setCoverageCrs( QString const &crs )
511 {
512 QgsDebugMsg( "Setting coverage CRS to " + crs + '.' );
513
514 if ( crs != mCoverageCrs && !crs.isEmpty() )
515 {
516 // delete old coordinate transform as it is no longer valid
517 mCoordinateTransform = QgsCoordinateTransform();
518
519 mExtentDirty = true;
520
521 mCoverageCrs = crs;
522
523 mCrs = QgsCoordinateReferenceSystem::fromOgcWmsCrs( mCoverageCrs );
524 }
525 }
526
setQueryItem(QUrl & url,const QString & item,const QString & value) const527 void QgsWcsProvider::setQueryItem( QUrl &url, const QString &item, const QString &value ) const
528 {
529 QUrlQuery query( url );
530 query.removeQueryItem( item );
531 query.addQueryItem( item, value );
532 url.setQuery( query );
533 }
534
readBlock(int bandNo,QgsRectangle const & viewExtent,int pixelWidth,int pixelHeight,void * block,QgsRasterBlockFeedback * feedback)535 bool QgsWcsProvider::readBlock( int bandNo, QgsRectangle const &viewExtent, int pixelWidth, int pixelHeight, void *block, QgsRasterBlockFeedback *feedback )
536 {
537 // TODO: set block to null values, move that to function and call only if fails
538 memset( block, 0, pixelWidth * pixelHeight * QgsRasterBlock::typeSize( dataType( bandNo ) ) );
539
540 // Requested extent must at least partially overlap coverage extent, otherwise
541 // server gives error. QGIS usually does not request blocks outside raster extent
542 // (higher level checks) but it is better to do check here as well
543 if ( !viewExtent.intersects( mCoverageExtent ) )
544 {
545 return false;
546 }
547
548 // Can we reuse the previously cached coverage?
549 if ( !mCachedGdalDataset ||
550 mCachedViewExtent != viewExtent ||
551 mCachedViewWidth != pixelWidth ||
552 mCachedViewHeight != pixelHeight )
553 {
554 getCache( bandNo, viewExtent, pixelWidth, pixelHeight, QString(), feedback );
555 }
556
557 if ( mCachedGdalDataset )
558 {
559 // It may happen (Geoserver) that if requested BBOX is larger than coverage
560 // extent, the returned data cover intersection of requested BBOX and coverage
561 // extent scaled to requested WIDTH/HEIGHT => check extent
562 // Unfortunately if the received raster does not have a CRS, the extent is the raster size
563 // and in that case it cannot be used to verify extent
564 QgsCoordinateReferenceSystem cacheCrs;
565 if ( !cacheCrs.createFromWkt( GDALGetProjectionRef( mCachedGdalDataset.get() ) ) &&
566 !cacheCrs.createFromWkt( GDALGetGCPProjection( mCachedGdalDataset.get() ) ) )
567 {
568 QgsDebugMsg( QStringLiteral( "Cached does not have CRS" ) );
569 }
570 QgsDebugMsg( "Cache CRS: " + cacheCrs.userFriendlyIdentifier() );
571
572 QgsRectangle cacheExtent = QgsGdalProviderBase::extent( mCachedGdalDataset.get() );
573 QgsDebugMsg( "viewExtent = " + viewExtent.toString() );
574 QgsDebugMsg( "cacheExtent = " + cacheExtent.toString() );
575 // TODO: check also rotated
576 if ( cacheCrs.isValid() && !mFixRotate )
577 {
578 // using qgsDoubleNear is too precise, example acceptable difference:
579 // 179.9999999306699863 x 179.9999999306700431
580 if ( !qgsDoubleNearSig( cacheExtent.xMinimum(), viewExtent.xMinimum(), 10 ) ||
581 !qgsDoubleNearSig( cacheExtent.yMinimum(), viewExtent.yMinimum(), 10 ) ||
582 !qgsDoubleNearSig( cacheExtent.xMaximum(), viewExtent.xMaximum(), 10 ) ||
583 !qgsDoubleNearSig( cacheExtent.yMaximum(), viewExtent.yMaximum(), 10 ) )
584 {
585 // Just print a message so user is aware of a server side issue but don't left
586 // the tile blank so we can deal with eventually misconfigured WCS server
587 // https://github.com/qgis/QGIS/issues/33339
588
589 QgsDebugMsg( QStringLiteral( "cacheExtent and viewExtent differ" ) );
590 QgsMessageLog::logMessage( tr( "Received coverage has wrong extent %1 (expected %2)" ).arg( cacheExtent.toString(), viewExtent.toString() ), tr( "WCS" ) );
591 }
592 }
593
594 int width = GDALGetRasterXSize( mCachedGdalDataset.get() );
595 int height = GDALGetRasterYSize( mCachedGdalDataset.get() );
596 QgsDebugMsg( QStringLiteral( "cached data width = %1 height = %2 (expected %3 x %4)" ).arg( width ).arg( height ).arg( pixelWidth ).arg( pixelHeight ) );
597
598 GDALRasterBandH gdalBand = GDALGetRasterBand( mCachedGdalDataset.get(), bandNo );
599 // TODO: check type?, check band count?
600 if ( mFixRotate && width == pixelHeight && height == pixelWidth )
601 {
602 // Rotate counter clockwise
603 // If GridOffsets With GeoServer,
604 QgsDebugMsg( QStringLiteral( "Rotating raster" ) );
605 int pixelSize = QgsRasterBlock::typeSize( dataType( bandNo ) );
606 QgsDebugMsg( QStringLiteral( "pixelSize = %1" ).arg( pixelSize ) );
607 int size = width * height * pixelSize;
608 void *tmpData = malloc( size );
609 if ( ! tmpData )
610 {
611 QgsDebugMsg( QStringLiteral( "Couldn't allocate memory of %1 bytes" ).arg( size ) );
612 return false;
613 }
614 if ( GDALRasterIO( gdalBand, GF_Read, 0, 0, width, height, tmpData, width, height, ( GDALDataType ) mGdalDataType.at( bandNo - 1 ), 0, 0 ) != CE_None )
615 {
616 QgsDebugMsg( QStringLiteral( "Raster IO Error" ) );
617 }
618 for ( int i = 0; i < pixelHeight; i++ )
619 {
620 for ( int j = 0; j < pixelWidth; j++ )
621 {
622 int destIndex = pixelSize * ( i * pixelWidth + j );
623 int srcIndex = pixelSize * ( j * width + ( width - i - 1 ) );
624 memcpy( ( char * )block + destIndex, ( char * )tmpData + srcIndex, pixelSize );
625 }
626 }
627 free( tmpData );
628 }
629 else if ( width == pixelWidth && height == pixelHeight )
630 {
631 if ( GDALRasterIO( gdalBand, GF_Read, 0, 0, pixelWidth, pixelHeight, block, pixelWidth, pixelHeight, ( GDALDataType ) mGdalDataType.at( bandNo - 1 ), 0, 0 ) != CE_None )
632 {
633 QgsDebugMsg( QStringLiteral( "Raster IO Error" ) );
634 }
635 else
636 {
637 QgsDebugMsg( QStringLiteral( "Block read OK" ) );
638 }
639 }
640 else
641 {
642 // This should not happen, but it is better to give distorted result + warning
643 if ( GDALRasterIO( gdalBand, GF_Read, 0, 0, width, height, block, pixelWidth, pixelHeight, ( GDALDataType ) mGdalDataType.at( bandNo - 1 ), 0, 0 ) != CE_None )
644 {
645 QgsDebugMsg( QStringLiteral( "Raster IO Error" ) );
646 }
647 QgsMessageLog::logMessage( tr( "Received coverage has wrong size %1 x %2 (expected %3 x %4)" ).arg( width ).arg( height ).arg( pixelWidth ).arg( pixelHeight ), tr( "WCS" ) );
648 }
649 }
650 return true;
651 }
652
getCache(int bandNo,QgsRectangle const & viewExtent,int pixelWidth,int pixelHeight,QString crs,QgsRasterBlockFeedback * feedback) const653 void QgsWcsProvider::getCache( int bandNo, QgsRectangle const &viewExtent, int pixelWidth, int pixelHeight, QString crs, QgsRasterBlockFeedback *feedback ) const
654 {
655 Q_UNUSED( bandNo )
656 // delete cached data
657 clearCache();
658
659 if ( crs.isEmpty() )
660 {
661 crs = mCoverageCrs;
662 }
663
664 mCachedViewExtent = viewExtent;
665 mCachedViewWidth = pixelWidth;
666 mCachedViewHeight = pixelHeight;
667
668 // --------------- BOUNDING BOX --------------------------------
669 //according to the WCS spec for 1.1, some CRS have inverted axis
670 // box:
671 // 1.0.0: minx,miny,maxx,maxy
672 // 1.1.0, 1.1.2: OGC 07-067r5 (WCS 1.1.2) refers to OGC 06-121r3 which says
673 // "The number of axes included, and the order of these axes, shall be as specified
674 // by the referenced CRS." That means inverted for geographic.
675 bool changeXY = false;
676 if ( !mIgnoreAxisOrientation && ( mCapabilities.version().startsWith( QLatin1String( "1.1" ) ) ) )
677 {
678 //create CRS from string
679 QgsCoordinateReferenceSystem srs = QgsCoordinateReferenceSystem::fromOgcWmsCrs( crs );
680 if ( srs.isValid() && srs.hasAxisInverted() )
681 {
682 changeXY = true;
683 }
684 }
685
686 if ( mInvertAxisOrientation ) changeXY = !changeXY;
687
688 double xRes = viewExtent.width() / pixelWidth;
689 double yRes = viewExtent.height() / pixelHeight;
690 QgsRectangle extent = viewExtent;
691 // WCS 1.1 grid is using grid points (surrounded by sample spaces) and
692 // "The spatial extent of a grid coverage extends only as far as the outermost
693 // grid points contained in the bounding box. It does NOT include any area
694 // (partial or whole grid cells or sample spaces) beyond those grid points."
695 // Mapserver and GDAL are using bbox defined by grid points, i.e. shrunk
696 // by 1 pixel, but Geoserver and ArcGIS are using full bbox including
697 // the space around edge grid points.
698 if ( mCapabilities.version().startsWith( QLatin1String( "1.1" ) ) && !mFixBox )
699 {
700 // shrink the extent to border cells centers by half cell size
701 extent = QgsRectangle( viewExtent.xMinimum() + xRes / 2., viewExtent.yMinimum() + yRes / 2., viewExtent.xMaximum() - xRes / 2., viewExtent.yMaximum() - yRes / 2. );
702 }
703
704 // Bounding box in WCS format (Warning: does not work with scientific notation)
705 QString bbox = QString( changeXY ? "%2,%1,%4,%3" : "%1,%2,%3,%4" )
706 .arg( qgsDoubleToString( extent.xMinimum() ),
707 qgsDoubleToString( extent.yMinimum() ),
708 qgsDoubleToString( extent.xMaximum() ),
709 qgsDoubleToString( extent.yMaximum() ) );
710
711 QUrl url( mIgnoreGetCoverageUrl ? mBaseUrl : mCapabilities.getCoverageUrl() );
712
713 // Version 1.0.0, 1.1.0, 1.1.2
714 setQueryItem( url, QStringLiteral( "SERVICE" ), QStringLiteral( "WCS" ) );
715 setQueryItem( url, QStringLiteral( "VERSION" ), mCapabilities.version() );
716 setQueryItem( url, QStringLiteral( "REQUEST" ), QStringLiteral( "GetCoverage" ) );
717 setQueryItem( url, QStringLiteral( "FORMAT" ), mFormat );
718
719 // Version 1.0.0
720 if ( mCapabilities.version().startsWith( QLatin1String( "1.0" ) ) )
721 {
722 setQueryItem( url, QStringLiteral( "COVERAGE" ), mIdentifier );
723 if ( !mTime.isEmpty() )
724 {
725 // It seems that Mmapserver (6.0.3) WCS 1.1 completely ignores
726 // TemporalDomain. Some code (copy-pasted from 1.0) is commented in
727 // msWCSDescribeCoverage_CoverageDescription11() and GetCoverage
728 // TimeSequence param is not supported at all. If a coverage is defined
729 // with timeposition in mapfile, the result of GetCoverage is empty
730 // raster (all values 0).
731 setQueryItem( url, QStringLiteral( "TIME" ), mTime );
732 }
733 setQueryItem( url, QStringLiteral( "BBOX" ), bbox );
734 setQueryItem( url, QStringLiteral( "CRS" ), crs ); // request BBOX CRS
735 setQueryItem( url, QStringLiteral( "RESPONSE_CRS" ), crs ); // response CRS
736 setQueryItem( url, QStringLiteral( "WIDTH" ), QString::number( pixelWidth ) );
737 setQueryItem( url, QStringLiteral( "HEIGHT" ), QString::number( pixelHeight ) );
738 }
739
740 // Version 1.1.0, 1.1.2
741 if ( mCapabilities.version().startsWith( QLatin1String( "1.1" ) ) )
742 {
743 setQueryItem( url, QStringLiteral( "IDENTIFIER" ), mIdentifier );
744 QString crsUrn = QStringLiteral( "urn:ogc:def:crs:%1::%2" ).arg( crs.split( ':' ).value( 0 ), crs.split( ':' ).value( 1 ) );
745 bbox += ',' + crsUrn;
746
747 if ( !mTime.isEmpty() )
748 {
749 setQueryItem( url, QStringLiteral( "TIMESEQUENCE" ), mTime );
750 }
751
752 setQueryItem( url, QStringLiteral( "BOUNDINGBOX" ), bbox );
753
754 // Example:
755 // GridBaseCRS=urn:ogc:def:crs:SG:6.6:32618
756 // GridType=urn:ogc:def:method:CS:1.1:2dGridIn2dCrs
757 // GridCS=urn:ogc:def:cs:OGC:0Grid2dSquareCS
758 // GridOrigin=0,0
759 // GridOffsets=0.0707,-0.0707,0.1414,0.1414&
760
761 setQueryItem( url, QStringLiteral( "GRIDBASECRS" ), crsUrn ); // response CRS
762
763 setQueryItem( url, QStringLiteral( "GRIDCS" ), QStringLiteral( "urn:ogc:def:cs:OGC:0.0:Grid2dSquareCS" ) );
764
765 setQueryItem( url, QStringLiteral( "GRIDTYPE" ), QStringLiteral( "urn:ogc:def:method:WCS:1.1:2dSimpleGrid" ) );
766
767 // GridOrigin is BBOX minx, maxy
768 // Note: shifting origin to cell center (not really necessary nor making sense)
769 // does not work with Mapserver 6.0.3
770 // Mapserver 6.0.3 does not work with origin on yMinimum (lower left)
771 // Geoserver works OK with yMinimum (lower left)
772 QString gridOrigin = QString( changeXY ? "%2,%1" : "%1,%2" )
773 .arg( qgsDoubleToString( extent.xMinimum() ),
774 qgsDoubleToString( extent.yMaximum() ) );
775 setQueryItem( url, QStringLiteral( "GRIDORIGIN" ), gridOrigin );
776
777 // GridOffsets WCS 1.1:
778 // GridType urn:ogc:def:method:WCS:1.1:2dSimpleGrid : 2 values
779 // GridType urn:ogc:def:method:WCS:1.1:2dGridIn2dCrs : 4 values, the center two of these offsets will be zero when the GridCRS is not rotated or skewed in the GridBaseCRS.
780 // The yOff must be negative because origin is in upper corner
781 // WCS 1.1.2: BaseX = origin(1) + offsets(1) * GridX
782 // BaseY = origin(2) + offsets(2) * GridY
783 // otherwise GeoServer gives mirrored result, however
784 // Mapserver works OK with both positive and negative
785 // OTOH, the yOff offset must not be negative with mFixRotate and Geoserver 2.1-SNAPSHOT
786 // but it must be negative with GeoServer 2.1.3 and mFixRotate. I am not sure
787 // at this moment 100% -> disabling positive yOff for now - TODO: try other servers
788 //double yOff = mFixRotate ? yRes : -yRes; // this was working with some servers I think
789 double yOff = -yRes;
790 QString gridOffsets = QString( changeXY ? "%2,%1" : "%1,%2" )
791 //setQueryItem( url, "GRIDTYPE", "urn:ogc:def:method:WCS:1.1:2dGridIn2dCrs" );
792 //QString gridOffsets = QString( changeXY ? "%2,0,0,%1" : "%1,0,0,%2" )
793 .arg( qgsDoubleToString( xRes ),
794 qgsDoubleToString( yOff ) );
795 setQueryItem( url, QStringLiteral( "GRIDOFFSETS" ), gridOffsets );
796 }
797
798 QgsDebugMsg( QStringLiteral( "GetCoverage: %1" ).arg( url.toString() ) );
799
800 // cache some details for if the user wants to do an identifyAsHtml() later
801
802 //mGetFeatureInfoUrlBase = mIgnoreGetFeatureInfoUrl ? mBaseUrl : getFeatureInfoUrl();
803
804 QgsWcsDownloadHandler handler( url, mAuth, mCacheLoadControl, mCachedData, mCapabilities.version(), mCachedError, feedback );
805 handler.blockingDownload();
806
807 QgsDebugMsg( QStringLiteral( "%1 bytes received" ).arg( mCachedData.size() ) );
808 if ( mCachedData.isEmpty() )
809 {
810 if ( !feedback || !feedback->isCanceled() )
811 {
812 QgsMessageLog::logMessage( tr( "No data received" ), tr( "WCS" ) );
813 }
814 clearCache();
815 return;
816 }
817
818 #if 0
819 QFile myFile( "/tmp/qgiswcscache.dat" );
820 if ( myFile.open( QIODevice::WriteOnly ) )
821 {
822 QDataStream myStream( &myFile );
823 myStream.writeRawData( mCachedData.data(), mCachedData.size() );
824 myFile.close();
825 }
826 #endif
827
828 mCachedMemFile = VSIFileFromMemBuffer( mCachedMemFilename.toUtf8().constData(),
829 ( GByte * )mCachedData.data(),
830 ( vsi_l_offset )mCachedData.size(),
831 FALSE );
832
833 if ( !mCachedMemFile )
834 {
835 QgsMessageLog::logMessage( tr( "Cannot create memory file" ), tr( "WCS" ) );
836 clearCache();
837 return;
838 }
839 QgsDebugMsg( QStringLiteral( "Memory file created" ) );
840
841 CPLErrorReset();
842 mCachedGdalDataset.reset( GDALOpen( mCachedMemFilename.toUtf8().constData(), GA_ReadOnly ) );
843 if ( !mCachedGdalDataset )
844 {
845 QgsMessageLog::logMessage( QString::fromUtf8( CPLGetLastErrorMsg() ), tr( "WCS" ) );
846 clearCache();
847 return;
848 }
849 QgsDebugMsg( QStringLiteral( "Dataset opened" ) );
850
851 }
852
853 // For stats only, maybe change QgsRasterDataProvider::bandStatistics() to
854 // use standard readBlock with extent
readBlock(int bandNo,int xBlock,int yBlock,void * block)855 bool QgsWcsProvider::readBlock( int bandNo, int xBlock, int yBlock, void *block )
856 {
857
858 QgsDebugMsg( QStringLiteral( "xBlock = %1 yBlock = %2" ).arg( xBlock ).arg( yBlock ) );
859
860 if ( !mHasSize )
861 return false;
862
863 double xRes = mCoverageExtent.width() / mWidth;
864 double yRes = mCoverageExtent.height() / mHeight;
865
866 // blocks on edges may run out of extent, that should not be problem (at least for
867 // stats - there is a check for it)
868 double xMin = mCoverageExtent.xMinimum() + xRes * xBlock * mXBlockSize;
869 double xMax = xMin + xRes * mXBlockSize;
870 double yMax = mCoverageExtent.yMaximum() - yRes * yBlock * mYBlockSize;
871 double yMin = yMax - yRes * mYBlockSize;
872 //QgsDebugMsg( QStringLiteral("yMin = %1 yMax = %2").arg(yMin).arg(yMax) );
873
874 QgsRectangle extent( xMin, yMin, xMax, yMax );
875
876 return readBlock( bandNo, extent, mXBlockSize, mYBlockSize, block, nullptr );
877 }
878
879
880 // This could be shared with GDAL provider
sourceDataType(int bandNo) const881 Qgis::DataType QgsWcsProvider::sourceDataType( int bandNo ) const
882 {
883 if ( bandNo <= 0 || bandNo > mSrcGdalDataType.size() )
884 {
885 return Qgis::UnknownDataType;
886 }
887
888 return dataTypeFromGdal( mSrcGdalDataType[bandNo - 1] );
889 }
890
dataType(int bandNo) const891 Qgis::DataType QgsWcsProvider::dataType( int bandNo ) const
892 {
893 if ( bandNo <= 0 || bandNo > mGdalDataType.size() )
894 {
895 return Qgis::UnknownDataType;
896 }
897
898 return dataTypeFromGdal( mGdalDataType[bandNo - 1] );
899 }
900
bandCount() const901 int QgsWcsProvider::bandCount() const
902 {
903 return mBandCount;
904 }
905
906 // this is only called once when statistics are calculated
907 // TODO
xBlockSize() const908 int QgsWcsProvider::xBlockSize() const
909 {
910 return mXBlockSize;
911 }
yBlockSize() const912 int QgsWcsProvider::yBlockSize() const
913 {
914 return mYBlockSize;
915 }
916
xSize() const917 int QgsWcsProvider::xSize() const { return mWidth; }
ySize() const918 int QgsWcsProvider::ySize() const { return mHeight; }
919
clearCache() const920 void QgsWcsProvider::clearCache() const
921 {
922 if ( mCachedGdalDataset )
923 {
924 QgsDebugMsg( QStringLiteral( "Close mCachedGdalDataset" ) );
925 mCachedGdalDataset.reset();
926 QgsDebugMsg( QStringLiteral( "Closed" ) );
927 }
928 if ( mCachedMemFile )
929 {
930 QgsDebugMsg( QStringLiteral( "Close mCachedMemFile" ) );
931 VSIFCloseL( mCachedMemFile );
932 mCachedMemFile = nullptr;
933 QgsDebugMsg( QStringLiteral( "Closed" ) );
934 }
935 QgsDebugMsg( QStringLiteral( "Clear mCachedData" ) );
936 mCachedData.clear();
937 mCachedError.clear();
938 QgsDebugMsg( QStringLiteral( "Cleared" ) );
939 }
940
colorTable(int bandNumber) const941 QList<QgsColorRampShader::ColorRampItem> QgsWcsProvider::colorTable( int bandNumber )const
942 {
943 return mColorTables.value( bandNumber - 1 );
944 }
945
colorInterpretation(int bandNo) const946 int QgsWcsProvider::colorInterpretation( int bandNo ) const
947 {
948 GDALRasterBandH myGdalBand = GDALGetRasterBand( mCachedGdalDataset.get(), bandNo );
949 return colorInterpretationFromGdal( GDALGetRasterColorInterpretation( myGdalBand ) );
950 }
951
952
parseServiceExceptionReportDom(const QByteArray & xml,const QString & wcsVersion,QString & errorTitle,QString & errorText)953 bool QgsWcsProvider::parseServiceExceptionReportDom( const QByteArray &xml, const QString &wcsVersion, QString &errorTitle, QString &errorText )
954 {
955
956 #ifdef QGISDEBUG
957 //test the content of the QByteArray
958 QString responsestring( xml );
959 QgsDebugMsg( "received the following data: " + responsestring );
960 #endif
961
962 // Convert completed document into a Dom
963 QDomDocument doc;
964 QString errorMsg;
965 int errorLine;
966 int errorColumn;
967 bool contentSuccess = doc.setContent( xml, false, &errorMsg, &errorLine, &errorColumn );
968
969 if ( !contentSuccess )
970 {
971 errorTitle = tr( "Dom Exception" );
972 errorText = tr( "Could not get WCS Service Exception at %1 at line %2 column %3\n\nResponse was:\n\n%4" )
973 .arg( errorMsg )
974 .arg( errorLine )
975 .arg( errorColumn )
976 .arg( QString( xml ) );
977
978 QgsLogger::debug( "Dom Exception: " + errorText );
979
980 return false;
981 }
982
983 QDomElement docElem = doc.documentElement();
984
985 // TODO: Assert the docElem.tagName() is
986 // ServiceExceptionReport // 1.0
987 // ows:ExceptionReport // 1.1
988
989 //QString version = docElem.attribute("version");
990
991 QDomElement e;
992 if ( wcsVersion.startsWith( QLatin1String( "1.0" ) ) )
993 {
994 e = QgsWcsCapabilities::domElement( docElem, QStringLiteral( "ServiceException" ) );
995 }
996 else // 1.1
997 {
998 e = QgsWcsCapabilities::domElement( docElem, QStringLiteral( "Exception" ) );
999 }
1000 parseServiceException( e, wcsVersion, errorTitle, errorText );
1001
1002 QgsDebugMsg( QStringLiteral( "exiting." ) );
1003
1004 return true;
1005 }
1006
parseServiceException(QDomElement const & e,const QString & wcsVersion,QString & errorTitle,QString & errorText)1007 void QgsWcsProvider::parseServiceException( QDomElement const &e, const QString &wcsVersion, QString &errorTitle, QString &errorText )
1008 {
1009
1010 errorTitle = tr( "Service Exception" );
1011
1012 QMap<QString, QString> exceptions;
1013
1014 // Some codes are in both 1.0 and 1.1, in that case the 'meaning of code' is
1015 // taken from 1.0 specification
1016
1017 // set up friendly descriptions for the service exception
1018 // 1.0
1019 exceptions[QStringLiteral( "InvalidFormat" )] = tr( "Request contains a format not offered by the server." );
1020 exceptions[QStringLiteral( "CoverageNotDefined" )] = tr( "Request is for a Coverage not offered by the service instance." );
1021 exceptions[QStringLiteral( "CurrentUpdateSequence" )] = tr( "Value of (optional) UpdateSequence parameter in GetCapabilities request is equal to current value of service metadata update sequence number." );
1022 exceptions[QStringLiteral( "InvalidUpdateSequence" )] = tr( "Value of (optional) UpdateSequence parameter in GetCapabilities request is greater than current value of service metadata update sequence number." );
1023 // 1.0, 1.1
1024 exceptions[QStringLiteral( "MissingParameterValue" )] = tr( "Request does not include a parameter value, and the server instance did not declare a default value for that dimension." );
1025 exceptions[QStringLiteral( "InvalidParameterValue" )] = tr( "Request contains an invalid parameter value." );
1026 // 1.1
1027 exceptions[QStringLiteral( "NoApplicableCode" )] = tr( "No other exceptionCode specified by this service and server applies to this exception." );
1028 exceptions[QStringLiteral( "UnsupportedCombination" )] = tr( "Operation request contains an output CRS that can not be used within the output format." );
1029 exceptions[QStringLiteral( "NotEnoughStorage" )] = tr( "Operation request specifies to \"store\" the result, but not enough storage is available to do this." );
1030
1031 QString seCode;
1032 QString seText;
1033 if ( wcsVersion.startsWith( QLatin1String( "1.0" ) ) )
1034 {
1035 seCode = e.attribute( QStringLiteral( "code" ) );
1036 seText = e.text();
1037 }
1038 else
1039 {
1040 QStringList codes;
1041 seCode = e.attribute( QStringLiteral( "exceptionCode" ) );
1042 // UMN Mapserver (6.0.3) has messed/switched 'locator' and 'exceptionCode'
1043 if ( ! exceptions.contains( seCode ) )
1044 {
1045 seCode = e.attribute( QStringLiteral( "locator" ) );
1046 if ( ! exceptions.contains( seCode ) )
1047 {
1048 seCode.clear();
1049 }
1050 }
1051 seText = QgsWcsCapabilities::firstChildText( e, QStringLiteral( "ExceptionText" ) );
1052 }
1053
1054 if ( seCode.isEmpty() )
1055 {
1056 errorText = tr( "(No error code was reported)" );
1057 }
1058 else if ( exceptions.contains( seCode ) )
1059 {
1060 errorText = exceptions.value( seCode );
1061 }
1062 else
1063 {
1064 errorText = seCode + ' ' + tr( "(Unknown error code)" );
1065 }
1066
1067 errorText += '\n' + tr( "The WCS vendor also reported: " );
1068 errorText += seText;
1069
1070 QgsMessageLog::logMessage( tr( "composed error message '%1'." ).arg( errorText ), tr( "WCS" ) );
1071 QgsDebugMsg( QStringLiteral( "exiting." ) );
1072 }
1073
1074
1075
extent() const1076 QgsRectangle QgsWcsProvider::extent() const
1077 {
1078 if ( mExtentDirty )
1079 {
1080 if ( calculateExtent() )
1081 {
1082 mExtentDirty = false;
1083 }
1084 }
1085
1086 return mCoverageExtent;
1087 }
1088
isValid() const1089 bool QgsWcsProvider::isValid() const
1090 {
1091 return mValid;
1092 }
1093
1094
wcsVersion()1095 QString QgsWcsProvider::wcsVersion()
1096 {
1097 return mCapabilities.version();
1098 }
1099
calculateExtent() const1100 bool QgsWcsProvider::calculateExtent() const
1101 {
1102
1103 // Make sure we know what extents are available
1104 if ( !mCoverageSummary.described )
1105 {
1106 return false;
1107 }
1108
1109 // Prefer to use extent from capabilities / coverage description because
1110 // transformation from WGS84 enlarges the extent
1111 mCoverageExtent = mCoverageSummary.boundingBoxes.value( mCoverageCrs );
1112 QgsDebugMsg( "mCoverageCrs = " + mCoverageCrs + " mCoverageExtent = " + mCoverageExtent.toString() );
1113 if ( !mCoverageExtent.isEmpty() && mCoverageExtent.isFinite() )
1114 {
1115 QgsDebugMsg( "mCoverageExtent = " + mCoverageExtent.toString() );
1116 //return true;
1117 }
1118 else
1119 {
1120 // Set up the coordinate transform from the WCS standard CRS:84 bounding
1121 // box to the user's selected CRS
1122 if ( !mCoordinateTransform.isValid() )
1123 {
1124 QgsCoordinateReferenceSystem qgisSrsSource = QgsCoordinateReferenceSystem::fromOgcWmsCrs( QStringLiteral( "EPSG:4326" ) );
1125 QgsCoordinateReferenceSystem qgisSrsDest = QgsCoordinateReferenceSystem::fromOgcWmsCrs( mCoverageCrs );
1126
1127 //QgsDebugMsg( "qgisSrsSource: " + qgisSrsSource.toWkt() );
1128 //QgsDebugMsg( "qgisSrsDest: " + qgisSrsDest.toWkt() );
1129
1130 mCoordinateTransform = QgsCoordinateTransform( qgisSrsSource, qgisSrsDest, transformContext() );
1131 }
1132
1133 QgsDebugMsg( "mCoverageSummary.wgs84BoundingBox= " + mCoverageSummary.wgs84BoundingBox.toString() );
1134
1135 // Convert to the user's CRS as required
1136 try
1137 {
1138 mCoverageExtent = mCoordinateTransform.transformBoundingBox( mCoverageSummary.wgs84BoundingBox, QgsCoordinateTransform::ForwardTransform );
1139 }
1140 catch ( QgsCsException &cse )
1141 {
1142 Q_UNUSED( cse )
1143 return false;
1144 }
1145
1146 //make sure extent does not contain 'inf' or 'nan'
1147 if ( !mCoverageExtent.isFinite() )
1148 {
1149 return false;
1150 }
1151 }
1152
1153 QgsDebugMsg( "mCoverageExtent = " + mCoverageExtent.toString() );
1154
1155 // It may happen (GeoServer) that extent reported in spatialDomain.Envelope is larger
1156 // than the coverage. Then if that larger BBOX is requested, the server returns
1157 // request BBOX intersected with coverage box scaled to requested WIDTH and HEIGHT.
1158 // GDAL WCS client does not suffer from this probably because it probably takes
1159 // extent from lonLatEnvelope (it probably does not calculate it from
1160 // spatialDomain.RectifiedGrid because calculated value is slightly different).
1161
1162 // To get the true extent (it can also be smaller than real if reported Envelope is
1163 // than real smaller, but smaller is safer because data cannot be shifted) we make
1164 // request of the whole extent cut the extent from spatialDomain.Envelope if
1165 // necessary
1166
1167 getCache( 1, mCoverageExtent, 10, 10 );
1168 if ( mCachedGdalDataset )
1169 {
1170 QgsRectangle cacheExtent = QgsGdalProviderBase::extent( mCachedGdalDataset.get() );
1171 QgsDebugMsg( "mCoverageExtent = " + mCoverageExtent.toString() );
1172 QgsDebugMsg( "cacheExtent = " + cacheExtent.toString() );
1173 QgsCoordinateReferenceSystem cacheCrs;
1174 if ( !cacheCrs.createFromWkt( GDALGetProjectionRef( mCachedGdalDataset.get() ) ) &&
1175 !cacheCrs.createFromWkt( GDALGetGCPProjection( mCachedGdalDataset.get() ) ) )
1176 {
1177 QgsDebugMsg( QStringLiteral( "Cached does not have CRS" ) );
1178 }
1179 QgsDebugMsg( "Cache CRS: " + cacheCrs.userFriendlyIdentifier() );
1180
1181 // We can only verify extent if CRS is set
1182 // If dataset comes rotated, GDAL probably cuts latitude extend, disable
1183 // extent check for rotated, TODO: verify
1184 if ( cacheCrs.isValid() && !mFixRotate )
1185 {
1186 if ( !qgsDoubleNearSig( cacheExtent.xMinimum(), mCoverageExtent.xMinimum(), 10 ) ||
1187 !qgsDoubleNearSig( cacheExtent.yMinimum(), mCoverageExtent.yMinimum(), 10 ) ||
1188 !qgsDoubleNearSig( cacheExtent.xMaximum(), mCoverageExtent.xMaximum(), 10 ) ||
1189 !qgsDoubleNearSig( cacheExtent.yMaximum(), mCoverageExtent.yMaximum(), 10 ) )
1190 {
1191 QgsDebugMsg( QStringLiteral( "cacheExtent and mCoverageExtent differ, mCoverageExtent cut to cacheExtent" ) );
1192 mCoverageExtent = cacheExtent;
1193 }
1194 }
1195 }
1196 else
1197 {
1198 // Unfortunately it may also happen that a server (cubewerx.com) does not have
1199 // overviews and it is not able to respond for the whole extent within timeout.
1200 // It returns timeout error.
1201 // In that case (if request failed) we do not report error to allow working
1202 // with such servers on smaller portions of extent
1203 // (http://lists.osgeo.org/pipermail/qgis-developer/2013-January/024019.html)
1204
1205 // Unfortunately even if we get over this 10x10 check, QGIS also requests
1206 // 32x32 thumbnail where it is waiting for another timeout
1207
1208 QgsDebugMsg( QStringLiteral( "Cannot get cache to verify extent" ) );
1209 QgsMessageLog::logMessage( tr( "Cannot verify coverage full extent: %1" ).arg( mCachedError.message() ), tr( "WCS" ) );
1210 }
1211
1212 return true;
1213 }
1214
1215
capabilities() const1216 int QgsWcsProvider::capabilities() const
1217 {
1218 int capability = NoCapabilities;
1219 capability |= QgsRasterDataProvider::Identify;
1220 capability |= QgsRasterDataProvider::IdentifyValue;
1221 capability |= QgsRasterDataProvider::Prefetch;
1222
1223 if ( mHasSize )
1224 {
1225 capability |= QgsRasterDataProvider::Size;
1226 }
1227
1228 return capability;
1229 }
1230
coverageMetadata(const QgsWcsCoverageSummary & coverage)1231 QString QgsWcsProvider::coverageMetadata( const QgsWcsCoverageSummary &coverage )
1232 {
1233 QString metadata;
1234
1235 // Use a nested table
1236 metadata += QLatin1String( "<tr><td>" );
1237 metadata += QLatin1String( "<table width=\"100%\">" );
1238
1239 // Table header
1240 metadata += QLatin1String( "<tr><th class=\"strong\">" );
1241 metadata += tr( "Property" );
1242 metadata += QLatin1String( "</th>" );
1243 metadata += QLatin1String( "<th class=\"strong\">" );
1244 metadata += tr( "Value" );
1245 metadata += QLatin1String( "</th></tr>" );
1246
1247 metadata += htmlRow( tr( "Name (identifier)" ), coverage.identifier );
1248 metadata += htmlRow( tr( "Title" ), coverage.title );
1249 metadata += htmlRow( tr( "Abstract" ), coverage.abstract );
1250
1251 if ( !coverage.metadataLink.metadataType.isNull() &&
1252 !coverage.metadataLink.xlinkHref.isNull() )
1253 {
1254 metadata += htmlRow( tr( "Metadata Type" ), coverage.metadataLink.metadataType );
1255 metadata += htmlRow( tr( "Metadata Link" ), coverage.metadataLink.xlinkHref );
1256 }
1257
1258 #if 0
1259 // We don't have size, nativeCrs, nativeBoundingBox etc. until describe coverage which would be heavy for all coverages
1260 metadata += htmlRow( tr( "Fixed Width" ), QString::number( coverage.width ) );
1261 metadata += htmlRow( tr( "Fixed Height" ), QString::number( coverage.height ) );
1262 metadata += htmlRow( tr( "Native CRS" ), coverage.nativeCrs );
1263 metadata += htmlRow( tr( "Native Bounding Box" ), coverage.nativeBoundingBox.toString() );
1264 #endif
1265
1266 metadata += htmlRow( tr( "WGS 84 Bounding Box" ), coverage.wgs84BoundingBox.toString() );
1267
1268 // Layer Coordinate Reference Systems
1269 // TODO(?): supportedCrs and supportedFormat are not available in 1.0
1270 // until coverage is described - it would be confusing to show it only if available
1271 #if 0
1272 for ( int j = 0; j < std::min( coverage.supportedCrs.size(), 10 ); j++ )
1273 {
1274 metadata += htmlRow( tr( "Available in CRS" ), coverage.supportedCrs.value( j ) );
1275 }
1276
1277 if ( coverage.supportedCrs.size() > 10 )
1278 {
1279 metadata += htmlRow( tr( "Available in CRS" ), tr( "(and %n more)", "crs", coverage.supportedCrs.size() - 10 ) );
1280 }
1281
1282 for ( int j = 0; j < std::min( coverage.supportedFormat.size(), 10 ); j++ )
1283 {
1284 metadata += htmlRow( tr( "Available in format" ), coverage.supportedFormat.value( j ) );
1285 }
1286
1287 if ( coverage.supportedFormat.size() > 10 )
1288 {
1289 metadata += htmlRow( tr( "Available in format" ), tr( "(and %n more)", "crs", coverage.supportedFormat.size() - 10 ) );
1290 }
1291 #endif
1292
1293 // Close the nested table
1294 metadata += QLatin1String( "</table>" );
1295 metadata += QLatin1String( "</td></tr>" );
1296
1297 return metadata;
1298 }
1299
htmlMetadata()1300 QString QgsWcsProvider::htmlMetadata()
1301 {
1302 QString metadata;
1303 metadata += QStringLiteral( "<tr><td class=\"highlight\">" ) + tr( "WCS Info" ) + QStringLiteral( "</td><td><div>" );
1304
1305 metadata += QLatin1String( "</a> <a href=\"#coverages\">" );
1306 metadata += tr( "Coverages" );
1307 metadata += QLatin1String( "</a>" );
1308
1309 #if 0
1310 // TODO
1311 metadata += "<a href=\"#cachestats\">";
1312 metadata += tr( "Cache Stats" );
1313 metadata += "</a> ";
1314 #endif
1315
1316 metadata += QLatin1String( "<br /><table class=\"tabular-view\">" ); // Nested table 1
1317
1318 // Server Properties section
1319 metadata += QLatin1String( "<tr><th class=\"strong\"><a name=\"serverproperties\"></a>" );
1320 metadata += tr( "Server Properties" );
1321 metadata += QLatin1String( "</th></tr>" );
1322
1323 // Use a nested table
1324 metadata += QLatin1String( "<tr><td>" );
1325 metadata += QLatin1String( "<table width=\"100%\">" ); // Nested table 2
1326
1327 // Table header
1328 metadata += QLatin1String( "<tr><th class=\"strong\">" );
1329 metadata += tr( "Property" );
1330 metadata += QLatin1String( "</th>" );
1331 metadata += QLatin1String( "<th class=\"strong\">" );
1332 metadata += tr( "Value" );
1333 metadata += QLatin1String( "</th></tr>" );
1334
1335 metadata += htmlRow( ( "WCS Version" ), mCapabilities.version() );
1336 metadata += htmlRow( tr( "Title" ), mCapabilities.capabilities().title );
1337 metadata += htmlRow( tr( "Abstract" ), mCapabilities.capabilities().abstract );
1338 #if 0
1339 // TODO: probably apply stylesheet in QgsWcsCapabilities and save as html
1340 metadata += htmlRow( tr( "Keywords" ), mCapabilities.service.keywordList.join( "<br />" ) );
1341 metadata += htmlRow( tr( "Online Resource" ), "-" );
1342 metadata += htmlRow( tr( "Contact Person" ),
1343 mCapabilities.service.contactInformation.contactPersonPrimary.contactPerson
1344 + "<br />" + mCapabilities.service.contactInformation.contactPosition;
1345 + "<br />" + mCapabilities.service.contactInformation.contactPersonPrimary.contactOrganization );
1346 metadata += htmlRow( tr( "Fees" ), mCapabilities.service.fees );
1347 metadata += htmlRow( tr( "Access Constraints" ), mCapabilities.service.accessConstraints );
1348 metadata += htmlRow( tr( "Image Formats" ), mCapabilities.capability.request.getMap.format.join( "<br />" ) );
1349 metadata += htmlRow( tr( "GetCapabilitiesUrl" ), mBaseUrl );
1350 #endif
1351 metadata += htmlRow( tr( "Get Coverage Url" ), mCapabilities.getCoverageUrl() + ( mIgnoreGetCoverageUrl ? tr( " <font color=\"red\">(advertised but ignored)</font>" ) : QString() ) );
1352
1353 // Close the nested table
1354 metadata += QLatin1String( "</table>" ); // End nested table 2
1355 metadata += QLatin1String( "</td></tr>" );
1356
1357 // Coverage properties
1358 metadata += QLatin1String( "<tr><th class=\"strong\"><a name=\"coverages\"></a>" );
1359 metadata += tr( "Coverages" );
1360 metadata += QLatin1String( "</th></tr>" );
1361
1362 // Dialog takes too long to open if there are too many coverages (1000 for example)
1363 int count = 0;
1364 const auto constCoverages = mCapabilities.coverages();
1365 for ( const QgsWcsCoverageSummary &c : constCoverages )
1366 {
1367 metadata += coverageMetadata( c );
1368 count++;
1369 if ( count >= 100 ) break;
1370 }
1371 metadata += QLatin1String( "</table>" );
1372 if ( count < mCapabilities.coverages().size() )
1373 {
1374 metadata += tr( "And %1 more coverages" ).arg( mCapabilities.coverages().size() - count );
1375 }
1376
1377 metadata += QLatin1String( "</table></div></td></tr>\n" ); // End nested table 1
1378 return metadata;
1379 }
1380
htmlCell(const QString & text)1381 QString QgsWcsProvider::htmlCell( const QString &text )
1382 {
1383 return "<td>" + text + "</td>";
1384 }
1385
htmlRow(const QString & text1,const QString & text2)1386 QString QgsWcsProvider:: htmlRow( const QString &text1, const QString &text2 )
1387 {
1388 return "<tr>" + htmlCell( text1 ) + htmlCell( text2 ) + "</tr>";
1389 }
1390
identify(const QgsPointXY & point,QgsRaster::IdentifyFormat format,const QgsRectangle & boundingBox,int width,int height,int)1391 QgsRasterIdentifyResult QgsWcsProvider::identify( const QgsPointXY &point, QgsRaster::IdentifyFormat format, const QgsRectangle &boundingBox, int width, int height, int /*dpi*/ )
1392 {
1393 QgsDebugMsg( QStringLiteral( "thePoint = %1 %2" ).arg( point.x(), 0, 'g', 10 ).arg( point.y(), 0, 'g', 10 ) );
1394 QgsDebugMsg( QStringLiteral( "theWidth = %1 height = %2" ).arg( width ).arg( height ) );
1395 QgsDebugMsg( "boundingBox = " + boundingBox.toString() );
1396 QMap<int, QVariant> results;
1397
1398 if ( format != QgsRaster::IdentifyFormatValue )
1399 {
1400 return QgsRasterIdentifyResult( QGS_ERROR( tr( "Format not supported" ) ) );
1401 }
1402
1403 if ( !extent().contains( point ) )
1404 {
1405 // Outside the raster
1406 for ( int i = 1; i <= bandCount(); i++ )
1407 {
1408 results.insert( i, QVariant() );
1409 }
1410 return QgsRasterIdentifyResult( QgsRaster::IdentifyFormatValue, results );
1411 }
1412
1413 QgsRectangle finalExtent = boundingBox;
1414 int maxSize = 2000;
1415 int cacheSize = 1000; // tile cache size if context is not defined or small
1416 double relResTol = 0.1; // relative resolution tolerance (10%)
1417
1418 // TODO: We are using cacheSize x cacheSize if context is not defined
1419 // or box is too small (in pixels). That is necessary to allow effective
1420 // map querying (valuetool for example). OTOH, it means reading of large
1421 // unnecessary amount of data for each identify() call if query points are too
1422 // distant (out of cache) for example when called from scripts
1423
1424 // if context size is to large we have to cut it, in that case caching big
1425 // big part does not make sense
1426 if ( finalExtent.isEmpty() || width == 0 || height == 0 ||
1427 width > maxSize || height > maxSize )
1428 {
1429 // context missing, use a small area around the point and highest resolution if known
1430
1431 double xRes, yRes;
1432 if ( mHasSize )
1433 {
1434 xRes = mCoverageExtent.width() / mWidth;
1435 yRes = mCoverageExtent.height() / mHeight;
1436 }
1437 else
1438 {
1439 // set resolution approximately to 1mm
1440 switch ( mCrs.mapUnits() )
1441 {
1442 case QgsUnitTypes::DistanceMeters:
1443 xRes = 0.001;
1444 break;
1445 case QgsUnitTypes::DistanceFeet:
1446 xRes = 0.003;
1447 break;
1448 case QgsUnitTypes::DistanceDegrees:
1449 // max length of degree of latitude on pole is 111694 m
1450 xRes = 1e-8;
1451 break;
1452 default:
1453 xRes = 0.001; // expecting meters
1454 }
1455 yRes = xRes;
1456 }
1457
1458 width = cacheSize;
1459 height = cacheSize;
1460
1461 finalExtent = QgsRectangle( point.x() - xRes * width / 2,
1462 point.y() - yRes * height / 2,
1463 point.x() + xRes * width / 2,
1464 point.y() + yRes * height / 2 );
1465
1466 double xResDiff = std::fabs( mCachedViewExtent.width() / mCachedViewWidth - xRes );
1467 double yResDiff = std::fabs( mCachedViewExtent.height() / mCachedViewHeight - yRes );
1468
1469 if ( !mCachedGdalDataset ||
1470 !mCachedViewExtent.contains( point ) ||
1471 mCachedViewWidth == 0 || mCachedViewHeight == 0 ||
1472 xResDiff / xRes > relResTol ||
1473 yResDiff / yRes > relResTol )
1474 {
1475 getCache( 1, finalExtent, width, height );
1476 }
1477 }
1478 else
1479 {
1480 // Use context -> effective caching (usually, if context is constant)
1481 QgsDebugMsg( QStringLiteral( "Using context extent and resolution" ) );
1482 // To use the cache it is sufficient to have point within cache and
1483 // similar resolution
1484 double xRes = finalExtent.width() / width;
1485 double yRes = finalExtent.height() / height;
1486 QgsDebugMsg( QStringLiteral( "width = %1 height = %2 xRes = %3 yRes = %4" ).arg( finalExtent.width() ).arg( finalExtent.height() ).arg( xRes ).arg( yRes ) );
1487 double xResDiff = std::fabs( mCachedViewExtent.width() / mCachedViewWidth - xRes );
1488 double yResDiff = std::fabs( mCachedViewExtent.height() / mCachedViewHeight - yRes );
1489 QgsDebugMsg( QStringLiteral( "xRes diff = %1 yRes diff = %2 relative xResDiff = %3 relative yResDiff = %4" ).arg( xResDiff ).arg( yResDiff ).arg( xResDiff / xRes ).arg( yResDiff / yRes ) );
1490 if ( !mCachedGdalDataset ||
1491 !mCachedViewExtent.contains( point ) ||
1492 xResDiff / xRes > relResTol ||
1493 yResDiff / yRes > relResTol )
1494 {
1495 // Identify map tool is now using fake context 1x1 pixel to get point/line
1496 // features from WMS. In such case we enlarge the extent to get data cached
1497 // for next queries around.
1498 // BTW: UMN Mapserver (6.0.3) seems to be buggy with 1x1 pixels request (returns 'no data' value
1499 if ( width < cacheSize )
1500 {
1501 int buffer = ( cacheSize - width ) / 2;
1502 width += 2 * buffer;
1503 finalExtent.setXMinimum( finalExtent.xMinimum() - xRes * buffer );
1504 finalExtent.setXMaximum( finalExtent.xMaximum() + xRes * buffer );
1505 }
1506 if ( height < cacheSize )
1507 {
1508 int buffer = ( cacheSize - height ) / 2;
1509 height += 2 * buffer;
1510 finalExtent.setYMinimum( finalExtent.yMinimum() - yRes * buffer );
1511 finalExtent.setYMaximum( finalExtent.yMaximum() + yRes * buffer );
1512 }
1513
1514 getCache( 1, finalExtent, width, height );
1515 }
1516 }
1517
1518 if ( !mCachedGdalDataset ||
1519 !mCachedViewExtent.contains( point ) )
1520 {
1521 return QgsRasterIdentifyResult( QGS_ERROR( tr( "Read data error" ) ) );
1522 }
1523
1524 double x = point.x();
1525 double y = point.y();
1526
1527 // Calculate the row / column where the point falls
1528 double xRes = mCachedViewExtent.width() / mCachedViewWidth;
1529 double yRes = mCachedViewExtent.height() / mCachedViewHeight;
1530
1531 // Offset, not the cell index -> flor
1532 int col = ( int ) std::floor( ( x - mCachedViewExtent.xMinimum() ) / xRes );
1533 int row = ( int ) std::floor( ( mCachedViewExtent.yMaximum() - y ) / yRes );
1534
1535 QgsDebugMsg( "row = " + QString::number( row ) + " col = " + QString::number( col ) );
1536
1537 for ( int i = 1; i <= GDALGetRasterCount( mCachedGdalDataset.get() ); i++ )
1538 {
1539 GDALRasterBandH gdalBand = GDALGetRasterBand( mCachedGdalDataset.get(), i );
1540
1541 double value;
1542 CPLErr err = GDALRasterIO( gdalBand, GF_Read, col, row, 1, 1,
1543 &value, 1, 1, GDT_Float64, 0, 0 );
1544
1545 if ( err != CPLE_None )
1546 {
1547 QgsLogger::warning( "RasterIO error: " + QString::fromUtf8( CPLGetLastErrorMsg() ) );
1548 return QgsRasterIdentifyResult( QGS_ERROR( tr( "RasterIO error: " ) + QString::fromUtf8( CPLGetLastErrorMsg() ) ) );
1549 }
1550
1551 // Apply no data and user no data
1552 if ( ( sourceHasNoDataValue( i ) && useSourceNoDataValue( i ) &&
1553 ( std::isnan( value ) || qgsDoubleNear( value, sourceNoDataValue( i ) ) ) ) ||
1554 ( QgsRasterRange::contains( value, userNoDataValues( i ) ) ) )
1555 {
1556 results.insert( i, QVariant() );
1557 }
1558 else
1559 {
1560 results.insert( i, value );
1561 }
1562 }
1563
1564 return QgsRasterIdentifyResult( QgsRaster::IdentifyFormatValue, results );
1565 }
1566
crs() const1567 QgsCoordinateReferenceSystem QgsWcsProvider::crs() const
1568 {
1569 return mCrs;
1570 }
1571
lastErrorTitle()1572 QString QgsWcsProvider::lastErrorTitle()
1573 {
1574 return mErrorCaption;
1575 }
1576
lastError()1577 QString QgsWcsProvider::lastError()
1578 {
1579 QgsDebugMsg( "returning '" + mError + "'." );
1580 return mError;
1581 }
1582
lastErrorFormat()1583 QString QgsWcsProvider::lastErrorFormat()
1584 {
1585 return mErrorFormat;
1586 }
1587
name() const1588 QString QgsWcsProvider::name() const
1589 {
1590 return WCS_KEY;
1591 }
1592
description() const1593 QString QgsWcsProvider::description() const
1594 {
1595 return WCS_DESCRIPTION;
1596 }
1597
reloadProviderData()1598 void QgsWcsProvider::reloadProviderData()
1599 {
1600 clearCache();
1601 }
1602
nodeAttribute(const QDomElement & e,const QString & name,const QString & defValue)1603 QString QgsWcsProvider::nodeAttribute( const QDomElement &e, const QString &name, const QString &defValue )
1604 {
1605 if ( e.hasAttribute( name ) )
1606 return e.attribute( name );
1607
1608 QDomNamedNodeMap map( e.attributes() );
1609 for ( int i = 0; i < map.size(); i++ )
1610 {
1611 QDomAttr attr( map.item( i ).toElement().toAttr() );
1612 if ( attr.name().compare( name, Qt::CaseInsensitive ) == 0 )
1613 return attr.value();
1614 }
1615
1616 return defValue;
1617 }
1618
supportedMimes()1619 QMap<QString, QString> QgsWcsProvider::supportedMimes()
1620 {
1621 QMap<QString, QString> mimes;
1622 GDALAllRegister();
1623
1624 QgsDebugMsg( QStringLiteral( "GDAL drivers cont %1" ).arg( GDALGetDriverCount() ) );
1625 for ( int i = 0; i < GDALGetDriverCount(); ++i )
1626 {
1627 GDALDriverH driver = GDALGetDriver( i );
1628 Q_CHECK_PTR( driver );
1629
1630 if ( !driver )
1631 {
1632 QgsLogger::warning( "unable to get driver " + QString::number( i ) );
1633 continue;
1634 }
1635
1636 QString desc = GDALGetDescription( driver );
1637
1638 QString mimeType = GDALGetMetadataItem( driver, "DMD_MIMETYPE", "" );
1639
1640 if ( mimeType.isEmpty() ) continue;
1641
1642 desc = desc.isEmpty() ? mimeType : desc;
1643
1644 QgsDebugMsg( "add GDAL format " + mimeType + ' ' + desc );
1645
1646 mimes[mimeType] = desc;
1647 }
1648 return mimes;
1649 }
1650
createProvider(const QString & uri,const QgsDataProvider::ProviderOptions & options,QgsDataProvider::ReadFlags flags)1651 QgsWcsProvider *QgsWcsProviderMetadata::createProvider( const QString &uri, const QgsDataProvider::ProviderOptions &options, QgsDataProvider::ReadFlags flags )
1652 {
1653 return new QgsWcsProvider( uri, options, flags );
1654 }
1655
1656 // ----------
1657
1658 int QgsWcsDownloadHandler::sErrors = 0;
1659
QgsWcsDownloadHandler(const QUrl & url,QgsWcsAuthorization & auth,QNetworkRequest::CacheLoadControl cacheLoadControl,QByteArray & cachedData,const QString & wcsVersion,QgsError & cachedError,QgsRasterBlockFeedback * feedback)1660 QgsWcsDownloadHandler::QgsWcsDownloadHandler( const QUrl &url, QgsWcsAuthorization &auth, QNetworkRequest::CacheLoadControl cacheLoadControl, QByteArray &cachedData, const QString &wcsVersion, QgsError &cachedError, QgsRasterBlockFeedback *feedback )
1661 : mAuth( auth )
1662 , mEventLoop( new QEventLoop )
1663 , mCachedData( cachedData )
1664 , mWcsVersion( wcsVersion )
1665 , mCachedError( cachedError )
1666 , mFeedback( feedback )
1667 {
1668 if ( feedback )
1669 {
1670 connect( feedback, &QgsFeedback::canceled, this, &QgsWcsDownloadHandler::canceled, Qt::QueuedConnection );
1671
1672 // rendering could have been canceled before we started to listen to canceled() signal
1673 // so let's check before doing the download and maybe quit prematurely
1674 if ( feedback->isCanceled() )
1675 return;
1676 }
1677
1678 QNetworkRequest request( url );
1679 QgsSetRequestInitiatorClass( request, QStringLiteral( "QgsWcsDownloadHandler" ) );
1680 if ( !mAuth.setAuthorization( request ) )
1681 {
1682 QgsMessageLog::logMessage( tr( "Network request update failed for authentication config" ),
1683 tr( "WCS" ) );
1684 return;
1685 }
1686 request.setAttribute( QNetworkRequest::CacheSaveControlAttribute, true );
1687 request.setAttribute( QNetworkRequest::CacheLoadControlAttribute, cacheLoadControl );
1688
1689 mCacheReply = QgsNetworkAccessManager::instance()->get( request );
1690 if ( !mAuth.setAuthorizationReply( mCacheReply ) )
1691 {
1692 mCacheReply->deleteLater();
1693 mCacheReply = nullptr;
1694 QgsMessageLog::logMessage( tr( "Network reply update failed for authentication config" ),
1695 tr( "WCS" ) );
1696 finish();
1697 return;
1698 }
1699 connect( mCacheReply, &QNetworkReply::finished, this, &QgsWcsDownloadHandler::cacheReplyFinished );
1700 connect( mCacheReply, &QNetworkReply::downloadProgress, this, &QgsWcsDownloadHandler::cacheReplyProgress );
1701 }
1702
~QgsWcsDownloadHandler()1703 QgsWcsDownloadHandler::~QgsWcsDownloadHandler()
1704 {
1705 delete mEventLoop;
1706 }
1707
blockingDownload()1708 void QgsWcsDownloadHandler::blockingDownload()
1709 {
1710 if ( mFeedback && mFeedback->isCanceled() )
1711 return; // nothing to do
1712
1713 mEventLoop->exec( QEventLoop::ExcludeUserInputEvents );
1714
1715 Q_ASSERT( !mCacheReply );
1716 }
1717
cacheReplyFinished()1718 void QgsWcsDownloadHandler::cacheReplyFinished()
1719 {
1720 QgsDebugMsg( QStringLiteral( "mCacheReply->error() = %1" ).arg( mCacheReply->error() ) );
1721 if ( mCacheReply->error() == QNetworkReply::NoError )
1722 {
1723 QVariant redirect = mCacheReply->attribute( QNetworkRequest::RedirectionTargetAttribute );
1724 if ( !redirect.isNull() )
1725 {
1726 mCacheReply->deleteLater();
1727
1728 QgsDebugMsg( QStringLiteral( "redirected getmap: %1" ).arg( redirect.toString() ) );
1729 QNetworkRequest request( redirect.toUrl() );
1730 QgsSetRequestInitiatorClass( request, QStringLiteral( "QgsWcsDownloadHandler" ) );
1731 if ( !mAuth.setAuthorization( request ) )
1732 {
1733 QgsMessageLog::logMessage( tr( "Network request update failed for authentication config" ),
1734 tr( "WCS" ) );
1735 return;
1736 }
1737 mCacheReply = QgsNetworkAccessManager::instance()->get( request );
1738 if ( !mAuth.setAuthorizationReply( mCacheReply ) )
1739 {
1740 mCacheReply->deleteLater();
1741 mCacheReply = nullptr;
1742 QgsMessageLog::logMessage( tr( "Network reply update failed for authentication config" ),
1743 tr( "WCS" ) );
1744 finish();
1745 return;
1746 }
1747 connect( mCacheReply, &QNetworkReply::finished, this, &QgsWcsDownloadHandler::cacheReplyFinished );
1748 connect( mCacheReply, &QNetworkReply::downloadProgress, this, &QgsWcsDownloadHandler::cacheReplyProgress );
1749
1750 return;
1751 }
1752
1753 QVariant status = mCacheReply->attribute( QNetworkRequest::HttpStatusCodeAttribute );
1754 QgsDebugMsg( QStringLiteral( "status = %1" ).arg( status.toInt() ) );
1755 if ( !status.isNull() && status.toInt() >= 400 )
1756 {
1757 QVariant phrase = mCacheReply->attribute( QNetworkRequest::HttpReasonPhraseAttribute );
1758
1759 QgsMessageLog::logMessage( tr( "Map request error (Status: %1; Reason phrase: %2; URL: %3)" )
1760 .arg( status.toInt() )
1761 .arg( phrase.toString(),
1762 mCacheReply->url().toString() ), tr( "WCS" ) );
1763
1764 mCacheReply->deleteLater();
1765 mCacheReply = nullptr;
1766
1767 finish();
1768 return;
1769 }
1770
1771 // Read response
1772
1773 QString contentType = mCacheReply->header( QNetworkRequest::ContentTypeHeader ).toString();
1774 QgsDebugMsg( "contentType: " + contentType );
1775
1776 // Exception
1777 // Content type examples: text/xml
1778 // application/vnd.ogc.se_xml;charset=UTF-8
1779 // application/xml
1780 if ( contentType.startsWith( QLatin1String( "text/" ), Qt::CaseInsensitive ) ||
1781 contentType.compare( QLatin1String( "application/xml" ), Qt::CaseInsensitive ) == 0 ||
1782 contentType.startsWith( QLatin1String( "application/vnd.ogc.se_xml" ), Qt::CaseInsensitive ) )
1783 {
1784 QString errorTitle, errorText;
1785 QByteArray text = mCacheReply->readAll();
1786 if ( ( contentType.compare( QLatin1String( "text/xml" ), Qt::CaseInsensitive ) == 0 ||
1787 contentType.compare( QLatin1String( "application/xml" ), Qt::CaseInsensitive ) == 0 ||
1788 contentType.startsWith( QLatin1String( "application/vnd.ogc.se_xml" ), Qt::CaseInsensitive ) )
1789 && QgsWcsProvider::parseServiceExceptionReportDom( text, mWcsVersion, errorTitle, errorText ) )
1790 {
1791 mCachedError.append( SRVERR( tr( "Map request error:<br>Title: %1<br>Error: %2<br>URL: <a href='%3'>%3</a>)" )
1792 .arg( errorTitle, errorText,
1793 mCacheReply->url().toString() ) ) );
1794 }
1795 else
1796 {
1797 QgsMessageLog::logMessage( tr( "Map request error (Status: %1; Response: %2; URL: %3)" )
1798 .arg( status.toInt() )
1799 .arg( QString::fromUtf8( text ),
1800 mCacheReply->url().toString() ), tr( "WCS" ) );
1801 }
1802
1803 mCacheReply->deleteLater();
1804 mCacheReply = nullptr;
1805
1806 finish();
1807 return;
1808 }
1809
1810 // WCS 1.1 gives response as multipart, e.g.
1811 if ( QgsNetworkReplyParser::isMultipart( mCacheReply ) )
1812 {
1813 QgsDebugMsg( QStringLiteral( "reply is multipart" ) );
1814 QgsNetworkReplyParser parser( mCacheReply );
1815
1816 if ( !parser.isValid() )
1817 {
1818 QgsMessageLog::logMessage( tr( "Cannot parse multipart response: %1" ).arg( parser.error() ), tr( "WCS" ) );
1819 mCacheReply->deleteLater();
1820 mCacheReply = nullptr;
1821
1822 finish();
1823 return;
1824 }
1825
1826 if ( parser.parts() < 2 )
1827 {
1828 QgsMessageLog::logMessage( tr( "Expected 2 parts, %1 received" ).arg( parser.parts() ), tr( "WCS" ) );
1829 mCacheReply->deleteLater();
1830 mCacheReply = nullptr;
1831
1832 finish();
1833 return;
1834 }
1835 else if ( parser.parts() > 2 )
1836 {
1837 // We will try the second one
1838 QgsMessageLog::logMessage( tr( "More than 2 parts (%1) received" ).arg( parser.parts() ), tr( "WCS" ) );
1839 }
1840
1841 QString transferEncoding = parser.rawHeader( 1, QStringLiteral( "Content-Transfer-Encoding" ).toLatin1() );
1842 QgsDebugMsg( "transferEncoding = " + transferEncoding );
1843
1844 // It may happen (GeoServer) that in part header is for example
1845 // Content-Type: image/tiff and Content-Transfer-Encoding: base64
1846 // but content is xml ExceptionReport which is not in base64
1847 QByteArray body = parser.body( 1 );
1848 if ( body.startsWith( "<?xml" ) )
1849 {
1850 QString errorTitle, errorText;
1851 if ( QgsWcsProvider::parseServiceExceptionReportDom( body, mWcsVersion, errorTitle, errorText ) )
1852 {
1853 QgsMessageLog::logMessage( tr( "Map request error (Title: %1; Error: %2; URL: %3)" )
1854 .arg( errorTitle, errorText,
1855 mCacheReply->url().toString() ), tr( "WCS" ) );
1856 }
1857 else
1858 {
1859 QgsMessageLog::logMessage( tr( "Map request error (Response: %1; URL: %2)" )
1860 .arg( QString::fromUtf8( body ),
1861 mCacheReply->url().toString() ), tr( "WCS" ) );
1862 }
1863
1864 mCacheReply->deleteLater();
1865 mCacheReply = nullptr;
1866
1867 finish();
1868 return;
1869 }
1870
1871 if ( transferEncoding == QLatin1String( "binary" ) )
1872 {
1873 mCachedData = body;
1874 }
1875 else if ( transferEncoding == QLatin1String( "base64" ) )
1876 {
1877 mCachedData = QByteArray::fromBase64( body );
1878 }
1879 else
1880 {
1881 QgsMessageLog::logMessage( tr( "Content-Transfer-Encoding %1 not supported" ).arg( transferEncoding ), tr( "WCS" ) );
1882 }
1883 }
1884 else
1885 {
1886 // Mime types for WCS 1.0 response should be usually image/<format>
1887 // (image/tiff, image/png, image/jpeg, image/png; mode=8bit, etc.)
1888 // but other mime types (like application/*) may probably also appear
1889
1890 // Create memory file
1891 mCachedData = mCacheReply->readAll();
1892 }
1893
1894 mCacheReply->deleteLater();
1895 mCacheReply = nullptr;
1896
1897 finish();
1898 }
1899 else
1900 {
1901 // report any errors except for the one we have caused by canceling the request
1902 if ( mCacheReply->error() != QNetworkReply::OperationCanceledError )
1903 {
1904 // Resend request if AlwaysCache
1905 QNetworkRequest request = mCacheReply->request();
1906 QgsSetRequestInitiatorClass( request, QStringLiteral( "QgsWcsDownloadHandler" ) );
1907 if ( request.attribute( QNetworkRequest::CacheLoadControlAttribute ).toInt() == QNetworkRequest::AlwaysCache )
1908 {
1909 QgsDebugMsg( QStringLiteral( "Resend request with PreferCache" ) );
1910 request.setAttribute( QNetworkRequest::CacheLoadControlAttribute, QNetworkRequest::PreferCache );
1911
1912 mCacheReply->deleteLater();
1913
1914 mCacheReply = QgsNetworkAccessManager::instance()->get( request );
1915 if ( !mAuth.setAuthorizationReply( mCacheReply ) )
1916 {
1917 mCacheReply->deleteLater();
1918 mCacheReply = nullptr;
1919 QgsMessageLog::logMessage( tr( "Network reply update failed for authentication config" ),
1920 tr( "WCS" ) );
1921 finish();
1922 return;
1923 }
1924 connect( mCacheReply, &QNetworkReply::finished, this, &QgsWcsDownloadHandler::cacheReplyFinished, Qt::DirectConnection );
1925 connect( mCacheReply, &QNetworkReply::downloadProgress, this, &QgsWcsDownloadHandler::cacheReplyProgress, Qt::DirectConnection );
1926
1927 return;
1928 }
1929
1930 sErrors++;
1931 if ( sErrors < 100 )
1932 {
1933 QgsMessageLog::logMessage( tr( "Map request failed [error: %1 url: %2]" ).arg( mCacheReply->errorString(), mCacheReply->url().toString() ), tr( "WCS" ) );
1934 }
1935 else if ( sErrors == 100 )
1936 {
1937 QgsMessageLog::logMessage( tr( "Not logging more than 100 request errors." ), tr( "WCS" ) );
1938 }
1939 }
1940
1941 mCacheReply->deleteLater();
1942 mCacheReply = nullptr;
1943
1944 finish();
1945 }
1946 }
1947
cacheReplyProgress(qint64 bytesReceived,qint64 bytesTotal)1948 void QgsWcsDownloadHandler::cacheReplyProgress( qint64 bytesReceived, qint64 bytesTotal )
1949 {
1950 Q_UNUSED( bytesReceived )
1951 Q_UNUSED( bytesTotal )
1952 QgsDebugMsgLevel( QStringLiteral( "%1 of %2 bytes of map downloaded." ).arg( bytesReceived ).arg( bytesTotal < 0 ? QString( "unknown number of" ) : QString::number( bytesTotal ) ), 3 );
1953 }
1954
canceled()1955 void QgsWcsDownloadHandler::canceled()
1956 {
1957 QgsDebugMsg( QStringLiteral( "Caught canceled() signal" ) );
1958 if ( mCacheReply )
1959 {
1960 QgsDebugMsg( QStringLiteral( "Aborting WCS network request" ) );
1961 mCacheReply->abort();
1962 }
1963 }
1964
1965
dataItemProviders() const1966 QList< QgsDataItemProvider * > QgsWcsProviderMetadata::dataItemProviders() const
1967 {
1968 QList<QgsDataItemProvider *> providers;
1969 providers << new QgsWcsDataItemProvider;
1970 return providers;
1971 }
1972
QgsWcsProviderMetadata()1973 QgsWcsProviderMetadata::QgsWcsProviderMetadata():
1974 QgsProviderMetadata( QgsWcsProvider::WCS_KEY, QgsWcsProvider::WCS_DESCRIPTION ) {}
1975
providerMetadataFactory()1976 QGISEXTERN QgsProviderMetadata *providerMetadataFactory()
1977 {
1978 return new QgsWcsProviderMetadata();
1979 }
1980