Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
identify reprojected raster: use 1x1 pixel context in layer crs to ge…
…t reasonable search tolerance for WMS GetFeatureInfo, fixes #5896
  • Loading branch information
blazek committed Jan 29, 2014
1 parent ddec663 commit 6207529
Show file tree
Hide file tree
Showing 3 changed files with 82 additions and 15 deletions.
21 changes: 17 additions & 4 deletions src/gui/qgsmaptoolidentify.cpp
Expand Up @@ -431,6 +431,7 @@ bool QgsMapToolIdentify::identifyRasterLayer( QList<IdentifyResult> *results, Qg
return false;
}

QgsPoint pointInCanvasCrs = point;
try
{
point = toLayerCoordinates( layer, point );
Expand Down Expand Up @@ -460,12 +461,24 @@ bool QgsMapToolIdentify::identifyRasterLayer( QList<IdentifyResult> *results, Qg
}

QgsRasterIdentifyResult identifyResult;
// We can only use context (extent, width, height) if layer is not reprojected,
// otherwise we don't know source resolution (size).
// We can only use current map canvas context (extent, width, height) if layer is not reprojected,
if ( mCanvas->hasCrsTransformEnabled() && dprovider->crs() != mCanvas->mapRenderer()->destinationCrs() )
{
viewExtent = toLayerCoordinates( layer, viewExtent );
identifyResult = dprovider->identify( point, format );
// To get some reasonable response for point/line WMS vector layers we must
// use a context with approximately a resolution in layer CRS units
// corresponding to current map canvas resolution (for examplei UMN Mapserver
// in msWMSFeatureInfo() -> msQueryByRect() is using requested pixel
// + TOLERANCE (layer param) for feature selection)
//
QgsRectangle r;
r.setXMinimum( pointInCanvasCrs.x() - mapUnitsPerPixel / 2. );
r.setXMaximum( pointInCanvasCrs.x() + mapUnitsPerPixel / 2. );
r.setYMinimum( pointInCanvasCrs.y() - mapUnitsPerPixel / 2. );
r.setYMaximum( pointInCanvasCrs.y() + mapUnitsPerPixel / 2. );
r = toLayerCoordinates( layer, r ); // will be a bit larger
// Mapserver (6.0.3, for example) does not work with 1x1 pixel box
// but that is fixed (the rect is enlarged) in the WMS provider
identifyResult = dprovider->identify( point, format, r, 1, 1 );
}
else
{
Expand Down
59 changes: 48 additions & 11 deletions src/providers/wcs/qgswcsprovider.cpp
Expand Up @@ -1619,6 +1619,15 @@ QgsRasterIdentifyResult QgsWcsProvider::identify( const QgsPoint & thePoint, Qgs

QgsRectangle myExtent = theExtent;
int maxSize = 2000;
int cacheSize = 1000; // tile cache size if context is not defined or small
double relResTol = 0.1; // relative resolution tolerance (10%)

// TODO: We are using cacheSize x cacheSize if context is not defined
// or box is too small (in pixels). That is necessary to allow effective
// map querying (valuetool for example). OTOH, it means reading of large
// unnecessary amount of data for each identify() call if query points are too
// distant (out of cache) for example when called from scripts

// if context size is to large we have to cut it, in that case caching big
// big part does not make sense
if ( myExtent.isEmpty() || theWidth == 0 || theHeight == 0 ||
Expand All @@ -1627,8 +1636,7 @@ QgsRasterIdentifyResult QgsWcsProvider::identify( const QgsPoint & thePoint, Qgs
// context missing, use a small area around the point and highest resolution if known

double xRes, yRes;
//if ( mHasSize )
if ( false )
if ( mHasSize )
{
xRes = mCoverageExtent.width() / mWidth;
yRes = mCoverageExtent.height() / mHeight;
Expand All @@ -1654,20 +1662,22 @@ QgsRasterIdentifyResult QgsWcsProvider::identify( const QgsPoint & thePoint, Qgs
yRes = xRes;
}

// 1000x1000 to cache data around
theWidth = 1000;
theHeight = 1000;
theWidth = cacheSize;
theHeight = cacheSize;

myExtent = QgsRectangle( thePoint.x() - xRes * theWidth / 2,
thePoint.y() - yRes * theHeight / 2,
thePoint.x() + xRes * theWidth / 2,
thePoint.y() + yRes * theHeight / 2 );

double xResDiff = qAbs( mCachedViewExtent.width() / mCachedViewWidth - xRes );
double yResDiff = qAbs( mCachedViewExtent.height() / mCachedViewHeight - yRes );

if ( !mCachedGdalDataset ||
!mCachedViewExtent.contains( thePoint ) ||
mCachedViewWidth == 0 || mCachedViewHeight == 0 ||
mCachedViewExtent.width() / mCachedViewWidth - xRes > TINY_VALUE ||
mCachedViewExtent.height() / mCachedViewHeight - yRes > TINY_VALUE )
xResDiff / xRes > relResTol ||
yResDiff / yRes > relResTol )
{
getCache( 1, myExtent, theWidth, theHeight );
}
Expand All @@ -1676,12 +1686,39 @@ QgsRasterIdentifyResult QgsWcsProvider::identify( const QgsPoint & thePoint, Qgs
{
// Use context -> effective caching (usually, if context is constant)
QgsDebugMsg( "Using context extent and resolution" );
// To use the cache it is sufficient to have point within cache and
// similar resolution
double xRes = myExtent.width() / theWidth;
double yRes = myExtent.height() / theHeight;
QgsDebugMsg( QString( "width = %1 height = %2 xRes = %3 yRes = %4" ).arg( myExtent.width() ).arg( myExtent.height() ).arg( xRes ).arg( yRes ) );
double xResDiff = qAbs( mCachedViewExtent.width() / mCachedViewWidth - xRes );
double yResDiff = qAbs( mCachedViewExtent.height() / mCachedViewHeight - yRes );
QgsDebugMsg( QString( "xRes diff = %1 yRes diff = %2 relative xResDiff = %3 relative yResDiff = %4" ).arg( xResDiff ).arg( yResDiff ).arg( xResDiff / xRes ).arg( yResDiff / yRes ) );
if ( !mCachedGdalDataset ||
mCachedViewExtent != theExtent ||
mCachedViewWidth != theWidth ||
mCachedViewHeight != theHeight )
!mCachedViewExtent.contains( thePoint ) ||
xResDiff / xRes > relResTol ||
yResDiff / yRes > relResTol )
{
getCache( 1, theExtent, theWidth, theHeight );
// Identify map tool is now using fake context 1x1 pixel to get point/line
// features from WMS. In such case we enlarge the extent to get data cached
// for next queries around.
// BTW: UMN Mapserver (6.0.3) seems to be buggy with 1x1 pixels request (returns 'no data' value
if ( theWidth < cacheSize )
{
int buffer = ( cacheSize - theWidth ) / 2;
theWidth += 2 * buffer;
myExtent.setXMinimum( myExtent.xMinimum() - xRes * buffer );
myExtent.setXMaximum( myExtent.xMaximum() + xRes * buffer );
}
if ( theHeight < cacheSize )
{
int buffer = ( cacheSize - theHeight ) / 2;
theHeight += 2 * buffer;
myExtent.setYMinimum( myExtent.yMinimum() - yRes * buffer );
myExtent.setYMaximum( myExtent.yMaximum() + yRes * buffer );
}

getCache( 1, myExtent, theWidth, theHeight );
}
}

Expand Down
17 changes: 17 additions & 0 deletions src/providers/wms/qgswmsprovider.cpp
Expand Up @@ -4167,6 +4167,9 @@ QgsRasterIdentifyResult QgsWmsProvider::identify( const QgsPoint & thePoint, Qgs
{
// We don't know original source resolution, so we take some small extent around the point.

// Warning: this does not work well with poin/line vector layers where search rectangle
// is based on pixel size (e.g. UMN Mapserver is using TOLERANCE layer param)

double xRes = 0.001; // expecting meters

// TODO: add CRS as class member
Expand Down Expand Up @@ -4208,6 +4211,20 @@ QgsRasterIdentifyResult QgsWmsProvider::identify( const QgsPoint & thePoint, Qgs
double xRes = myExtent.width() / theWidth;
double yRes = myExtent.height() / theHeight;

// Mapserver (6.0.3, for example) does not seem to work with 1x1 pixel box
// (seems to be a different issue, not the slownes of GDAL with ECW mentioned above)
// so we have to enlarge it a bit
if ( theWidth == 1 )
{
theWidth += 1;
myExtent.setXMaximum( myExtent.xMaximum() + xRes );
}
if ( theHeight == 1 )
{
theHeight += 1;
myExtent.setYMaximum( myExtent.yMaximum() + yRes );
}

QgsDebugMsg( "myExtent = " + myExtent.toString() );
QgsDebugMsg( QString( "theWidth = %1 theHeight = %2" ).arg( theWidth ).arg( theHeight ) );
QgsDebugMsg( QString( "xRes = %1 yRes = %2" ).arg( xRes ).arg( yRes ) );
Expand Down

0 comments on commit 6207529

Please sign in to comment.