Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
wcs extent verification
  • Loading branch information
blazek committed Sep 25, 2012
1 parent f6e1fd3 commit a24f9c5
Show file tree
Hide file tree
Showing 4 changed files with 129 additions and 31 deletions.
29 changes: 29 additions & 0 deletions src/providers/gdal/qgsgdalproviderbase.cpp
Expand Up @@ -247,3 +247,32 @@ void QgsGdalProviderBase::registerGdalDrivers()
QgsApplication::applyGdalSkippedDrivers();
}
}

QgsRectangle QgsGdalProviderBase::extent( GDALDatasetH gdalDataset )const
{
double myGeoTransform[6];

bool myHasGeoTransform = GDALGetGeoTransform( gdalDataset, myGeoTransform ) == CE_None;
if ( !myHasGeoTransform )
{
// Initialise the affine transform matrix
myGeoTransform[0] = 0;
myGeoTransform[1] = 1;
myGeoTransform[2] = 0;
myGeoTransform[3] = 0;
myGeoTransform[4] = 0;
myGeoTransform[5] = -1;
}

// Use the affine transform to get geo coordinates for
// the corners of the raster
double myXMax = myGeoTransform[0] +
GDALGetRasterXSize( gdalDataset ) * myGeoTransform[1] +
GDALGetRasterYSize( gdalDataset ) * myGeoTransform[2];
double myYMin = myGeoTransform[3] +
GDALGetRasterXSize( gdalDataset ) * myGeoTransform[4] +
GDALGetRasterYSize( gdalDataset ) * myGeoTransform[5];

QgsRectangle extent( myGeoTransform[0], myYMin, myXMax, myGeoTransform[3] );
return extent;
}
2 changes: 2 additions & 0 deletions src/providers/gdal/qgsgdalproviderbase.h
Expand Up @@ -53,6 +53,8 @@ class QgsGdalProviderBase
int colorInterpretationFromGdal( int gdalColorInterpretation ) const;

QList<QgsColorRampShader::ColorRampItem> colorTable( GDALDatasetH gdalDataset, int bandNo )const;

QgsRectangle extent( GDALDatasetH gdalDataset ) const;
};

#endif
Expand Down
114 changes: 84 additions & 30 deletions src/providers/wcs/qgswcsprovider.cpp
Expand Up @@ -97,6 +97,7 @@ QgsWcsProvider::QgsWcsProvider( QString const &uri )
QgsDebugMsg( "constructing with uri '" + mHttpUri + "'." );

mValid = false;
mCachedMemFilename = QString( "/vsimem/qgis/wcs/%0.dat" ).arg(( qlonglong )this );

if ( !parseUri( uri ) ) return;

Expand Down Expand Up @@ -186,8 +187,6 @@ QgsWcsProvider::QgsWcsProvider( QString const &uri )
return;
}

mCachedMemFilename = QString( "/vsimem/qgis/wcs/%0.dat" ).arg(( qlonglong )this );

// Get small piece of coverage to find GDAL data type and nubmer of bands
int bandNo = 0; // All bands
int width;
Expand Down Expand Up @@ -526,9 +525,29 @@ void QgsWcsProvider::readBlock( int bandNo, QgsRectangle const & viewExtent, in

if ( mCachedGdalDataset )
{
// TODO band
int width = GDALGetRasterXSize( mCachedGdalDataset );
int height = GDALGetRasterYSize( mCachedGdalDataset );

// It may happen (Geoserver) that if requested BBOX is larger than coverage
// extent, the returned data cover intersection of requested BBOX and coverage
// extent scaled to requested WIDTH/HEIGHT => check extent
QgsRectangle cacheExtent = QgsGdalProviderBase::extent( mCachedGdalDataset );
QgsDebugMsg( "viewExtent = " + viewExtent.toString() );
QgsDebugMsg( "cacheExtent = " + cacheExtent.toString() );
if ( !doubleNear( cacheExtent.xMinimum(), viewExtent.xMinimum() ) ||
!doubleNear( cacheExtent.yMinimum(), viewExtent.yMinimum() ) ||
!doubleNear( cacheExtent.xMaximum(), viewExtent.xMaximum() ) ||
!doubleNear( cacheExtent.yMaximum(), viewExtent.yMaximum() ) )
{
QgsDebugMsg( "cacheExtent and viewExtent differ" );
QgsMessageLog::logMessage( tr( "Received coverage has wrong extent %1 (expected %2)" ).arg( cacheExtent.toString() ).arg( viewExtent.toString() ), tr( "WCS" ) );
// We are doing all possible to avoid this situation,
// If it happens, it would be possible to rescale the portion we get
// to only part of the data block, but it is better to left it
// blank, so that the problem may be discovered in its origin.
return;
}

GDALRasterBandH gdalBand = GDALGetRasterBand( mCachedGdalDataset, bandNo );
QgsDebugMsg( QString( "cached data width = %1 height = %2 (expected %3 x %4)" ).arg( width ).arg( height ).arg( pixelWidth ).arg( pixelHeight ) );
// TODO: check type? , check band count?
Expand Down Expand Up @@ -1295,52 +1314,87 @@ bool QgsWcsProvider::calculateExtent()
}

// Prefer to use extent from capabilities / coverage description because
// transformation from WGS84 increases the extent
// transformation from WGS84 enlarges the extent
mCoverageExtent = mCoverageSummary.boundingBoxes.value( mCoverageCrs );
QgsDebugMsg( "mCoverageCrs = " + mCoverageCrs + " mCoverageExtent = " + mCoverageExtent.toString() );
if ( !mCoverageExtent.isEmpty() && mCoverageExtent.isFinite() )
{
QgsDebugMsg( "mCoverageExtent = " + mCoverageExtent.toString() );
return true;
//return true;
}

// Set up the coordinate transform from the WCS standard CRS:84 bounding
// box to the user's selected CRS
if ( !mCoordinateTransform )
else
{
QgsCoordinateReferenceSystem qgisSrsSource;
QgsCoordinateReferenceSystem qgisSrsDest;
// Set up the coordinate transform from the WCS standard CRS:84 bounding
// box to the user's selected CRS
if ( !mCoordinateTransform )
{
QgsCoordinateReferenceSystem qgisSrsSource;
QgsCoordinateReferenceSystem qgisSrsDest;

//qgisSrsSource.createFromOgcWmsCrs( DEFAULT_LATLON_CRS );
qgisSrsSource.createFromOgcWmsCrs( "EPSG:4326" );
//QgsDebugMsg( "qgisSrsSource: " + qgisSrsSource.toWkt() );
qgisSrsDest.createFromOgcWmsCrs( mCoverageCrs );
//QgsDebugMsg( "qgisSrsDest: " + qgisSrsDest.toWkt() );
//qgisSrsSource.createFromOgcWmsCrs( DEFAULT_LATLON_CRS );
qgisSrsSource.createFromOgcWmsCrs( "EPSG:4326" );
//QgsDebugMsg( "qgisSrsSource: " + qgisSrsSource.toWkt() );
qgisSrsDest.createFromOgcWmsCrs( mCoverageCrs );
//QgsDebugMsg( "qgisSrsDest: " + qgisSrsDest.toWkt() );

mCoordinateTransform = new QgsCoordinateTransform( qgisSrsSource, qgisSrsDest );
mCoordinateTransform = new QgsCoordinateTransform( qgisSrsSource, qgisSrsDest );

}
}

QgsDebugMsg( "mCoverageSummary.wgs84BoundingBox= " + mCoverageSummary.wgs84BoundingBox.toString() );
QgsDebugMsg( "mCoverageSummary.wgs84BoundingBox= " + mCoverageSummary.wgs84BoundingBox.toString() );

// Convert to the user's CRS as required
try
{
mCoverageExtent = mCoordinateTransform->transformBoundingBox( mCoverageSummary.wgs84BoundingBox, QgsCoordinateTransform::ForwardTransform );
// Convert to the user's CRS as required
try
{
mCoverageExtent = mCoordinateTransform->transformBoundingBox( mCoverageSummary.wgs84BoundingBox, QgsCoordinateTransform::ForwardTransform );
}
catch ( QgsCsException &cse )
{
Q_UNUSED( cse );
return false;
}

//make sure extent does not contain 'inf' or 'nan'
if ( !mCoverageExtent.isFinite() )
{
return false;
}
}
catch ( QgsCsException &cse )

QgsDebugMsg( "mCoverageExtent = " + mCoverageExtent.toString() );

// It may happen (GeoServer) that extent reported in spatialDomain.Envelope is larger // than the coverage. Then if that larger BBOX is requested, the server returns
// request BBOX intersected with coverage box scaled to requested WIDTH and HEIGHT.
// GDAL WCS client does not suffer from this probably because it probably takes
// extent from lonLatEnvelope (it probably does not calculate it from
// spatialDomain.RectifiedGrid because calculated value is slightly different).

// To get the true extent (it can also be smaller than real if reported Envelope is
// than real smaller, but smaller is safer because data cannot be shifted) we make
// request of the whole extent cut the extent from spatialDomain.Envelope if
// necessary

getCache( 1, mCoverageExtent, 10, 10 );
if ( mCachedGdalDataset )
{
Q_UNUSED( cse );
return false;
QgsRectangle cacheExtent = QgsGdalProviderBase::extent( mCachedGdalDataset );
QgsDebugMsg( "mCoverageExtent = " + mCoverageExtent.toString() );
QgsDebugMsg( "cacheExtent = " + cacheExtent.toString() );
if ( !doubleNear( cacheExtent.xMinimum(), mCoverageExtent.xMinimum() ) ||
!doubleNear( cacheExtent.yMinimum(), mCoverageExtent.yMinimum() ) ||
!doubleNear( cacheExtent.xMaximum(), mCoverageExtent.xMaximum() ) ||
!doubleNear( cacheExtent.yMaximum(), mCoverageExtent.yMaximum() ) )
{
QgsDebugMsg( "cacheExtent and mCoverageExtent differ, mCoverageExtent cut to cacheExtent" );
mCoverageExtent = cacheExtent;
}
}

//make sure extent does not contain 'inf' or 'nan'
if ( !mCoverageExtent.isFinite() )
else
{
QgsDebugMsg( "Cannot get cache to verify extent" );
return false;
}

QgsDebugMsg( "mCoverageExtent = " + mCoverageExtent.toString() );
return true;
}

Expand Down
15 changes: 14 additions & 1 deletion tests/src/providers/wcs-servers.json
Expand Up @@ -41,7 +41,20 @@
}
]
}, {
url: 'http://demo.geonode.org/geoserver/wcs'
url: 'http://demo.geonode.org/geoserver/wcs',
issues: [
{
offender: 'server',
coverages: [ 'geonode:__137_FMCvivo_tif', 'geonode:__137_FMCvivo_tif_1', 'geonode:cl_minu1_2', 'geonode:cl_minu1_2_1', 'geonode:cl_minu1_2_2', 'geonode:dof_5modelAvg_sresb1_2010_2019', 'geonode:dof_5modelAvg_sresb1_2010_2019_1', 'geonode:MOS_EU_LAEA_2000' ],
versions: [ '1.0.0' ],
description: "The server DescribeCoverage fails: 'java.io.IOException null Translator error Unexpected error occurred during describe coverage xml encodin...'"
},{
offender: 'server',
coverages: [ 'geonode:elevation_1', 'geonode:elevation_2', 'geonode:quake' ],
versions: [ '1.0.0' ],
description: "The spatialDomain.Envelope in DescribeCoverage is wrong, it is larger than (real?) extent and that reported in lonLatEnvelope."
}
]
}, {
url: 'http://demo.mapserver.org/cgi-bin/wcs',
issues: [
Expand Down

0 comments on commit a24f9c5

Please sign in to comment.