WMS identify fixes

Author: blazek

src/app/qgsmaptoolidentify.cpp

@@ -379,17 +379,29 @@ bool QgsMapToolIdentify::identifyRasterLayer( QgsRasterLayer *layer, int x, int
     viewExtent = toLayerCoordinates( layer, viewExtent );
   }
 
-  // cut by layer extent to use possible cache, the same done when drawing
-  viewExtent = dprovider->extent().intersect( &viewExtent );
+  // It would be nice to use the same extent and size which was used for drawing,
+  // so that WCS can use cache from last draw, unfortunately QgsRasterLayer::draw()
+  // is doing some tricks with extent and size to allign raster to output which
+  // would be difficult to replicate here.
+  // Note: cutting the extent may result in slightly different x and y resolutions
+  // and thus shifted point calculated back in QGIS WMS (using average resolution)
+  //viewExtent = dprovider->extent().intersect( &viewExtent );
 
   double mapUnitsPerPixel = mCanvas->mapUnitsPerPixel();
   // Width and height are calculated from not projected extent and we hope that
   // are similar to source width and height used to reproject layer for drawing.
-  int width = mCanvas->extent().width() / mapUnitsPerPixel;
-  int height = mCanvas->extent().height() / mapUnitsPerPixel;
+  // TODO: may be very dangerous, because it may result in different resolutions
+  // in source CRS, and WMS server (QGIS server) calcs wrong coor using average resolution.
+  int width = qRound( viewExtent.width() / mapUnitsPerPixel );
+  int height = qRound( viewExtent.height() / mapUnitsPerPixel );
+
+  QgsDebugMsg( QString( "viewExtent.width = %1 viewExtent.height = %2" ).arg( viewExtent.width() ).arg( viewExtent.height() ) );
+  QgsDebugMsg( QString( "width = %1 height = %2" ).arg( width ).arg( height ) );
+  QgsDebugMsg( QString( "xRes = %1 yRes = %2 mapUnitsPerPixel = %3" ).arg( viewExtent.width() / width ).arg( viewExtent.height() / height ).arg( mapUnitsPerPixel ) );
 
   QMap< QString, QString > attributes, derivedAttributes;
 
+  QgsDebugMsg( QString( "idPoint: %1 %2" ).arg( idPoint.x() ).arg( idPoint.y() ) );
   attributes = dprovider->identify( idPoint, viewExtent, width, height );
 
   QString type;

src/core/qgsmaprenderer.cpp

@@ -206,6 +206,10 @@ void QgsMapRenderer::adjustExtentToSize()
   mExtent.setYMinimum( dymin );
   mExtent.setYMaximum( dymax );
 
+  QgsDebugMsg( QString( "Adjusted map units per pixel (x,y) : %1, %2" ).arg( mExtent.width() / myWidth, 0, 'f', 8 ).arg( mExtent.height() / myHeight, 0, 'f', 8 ) );
+
+  QgsDebugMsg( QString( "Recalced pixmap dimensions (x,y) : %1, %2" ).arg( mExtent.width() / mMapUnitsPerPixel, 0, 'f', 8 ).arg( mExtent.height() / mMapUnitsPerPixel, 0, 'f', 8 ) );
+
   // update the scale
   updateScale();
 

src/mapserver/qgspostrequesthandler.cpp

@@ -14,7 +14,7 @@
  *   (at your option) any later version.                                   *
  *                                                                         *
  ***************************************************************************/
-
+#include <stdlib.h>
 #include "qgspostrequesthandler.h"
 #include "qgslogger.h"
 #include <QDomDocument>
@@ -54,7 +54,7 @@ QMap<QString, QString> QgsPostRequestHandler::parseInput()
       QgsDebugMsg( "error, no query string found but a QDomDocument" );
       return parameters; //no query string? something must be wrong...
     }
-    
+
     requestStringToParameterMap( queryString, parameters );
 
     QDomElement docElem = doc.documentElement();

src/mapserver/qgswmsserver.cpp

@@ -648,6 +648,10 @@ int QgsWMSServer::getFeatureInfo( QDomDocument& result, QString version )
     return 2;
   }
 
+  QgsDebugMsg( "mMapRenderer->extent(): " +  mMapRenderer->extent().toString() );
+  QgsDebugMsg( QString( "mMapRenderer width = %1 height = %2" ).arg( mMapRenderer->outputSize().width() ).arg( mMapRenderer->outputSize().height() ) );
+  QgsDebugMsg( QString( "mMapRenderer->mapUnitsPerPixel() = %1" ).arg( mMapRenderer->mapUnitsPerPixel() ) );
+
   //find out the current scale denominater and set it to the SLD parser
   QgsScaleCalculator scaleCalc(( outputImage->logicalDpiX() + outputImage->logicalDpiY() ) / 2 , mMapRenderer->destinationCrs().mapUnits() );
   QgsRectangle mapExtent = mMapRenderer->extent();
@@ -1127,9 +1131,28 @@ int QgsWMSServer::infoPointToLayerCoordinates( int i, int j, QgsPoint* layerCoor
   }
 
   //first transform i,j to map output coordinates
-  QgsPoint mapPoint = mapRender->coordinateTransform()->toMapCoordinates( i, j );
+  // toMapCoordinates() is currently (Oct 18 2012) using average resolution
+  // to calc point but GetFeatureInfo request may be sent with different
+  // resolutions in each axis
+  //QgsPoint mapPoint = mapRender->coordinateTransform()->toMapCoordinates( i, j );
+  double xRes = mapRender->extent().width() / mapRender->width();
+  double yRes = mapRender->extent().height() / mapRender->height();
+  QgsPoint mapPoint( mapRender->extent().xMinimum() + i * xRes,
+                     mapRender->extent().yMaximum() - j * yRes );
+
+  QgsDebugMsg( QString( "mapPoint (corner): %1 %2" ).arg( mapPoint.x() ).arg( mapPoint.y() ) );
+  // use pixel center instead of corner
+  // Unfortunately going through pixel (integer) we cannot reconstruct precisely
+  // the coordinate clicked on client and thus result may differ from
+  // the same raster loaded and queried localy on client
+  mapPoint.setX( mapPoint.x() + xRes / 2 );
+  mapPoint.setY( mapPoint.y() - yRes / 2 );
+
+  QgsDebugMsg( QString( "mapPoint (pixel center): %1 %2" ).arg( mapPoint.x() ).arg( mapPoint.y() ) );
+
   //and then to layer coordinates
   *layerCoords = mapRender->mapToLayerCoordinates( layer, mapPoint );
+  QgsDebugMsg( QString( "mapPoint: %1 %2" ).arg( mapPoint.x() ).arg( mapPoint.y() ) );
   return 0;
 }
 
@@ -1303,9 +1326,11 @@ int QgsWMSServer::featureInfoFromRasterLayer( QgsRasterLayer* layer,
     return 1;
   }
 
+  QgsDebugMsg( QString( "infoPoint: %1 %2" ).arg( infoPoint->x() ).arg( infoPoint->y() ) );
+
   QMap<QString, QString> attributes;
-  // TODO: use context extent, width height (comes with request) to use WCS cache
-  attributes = layer->dataProvider()->identify( *infoPoint );
+  // use context extent, width height (comes with request) to use WCS cache
+  attributes = layer->dataProvider()->identify( *infoPoint, mMapRenderer->extent(), mMapRenderer->outputSize().width(), mMapRenderer->outputSize().height() );
 
   for ( QMap<QString, QString>::const_iterator it = attributes.constBegin(); it != attributes.constEnd(); ++it )
   {

src/providers/gdal/qgsgdalprovider.cpp

@@ -834,7 +834,7 @@ QMap<int, QVariant> QgsGdalProvider::identify( const QgsPoint & thePoint, Identi
   Q_UNUSED( theExtent );
   Q_UNUSED( theWidth );
   Q_UNUSED( theHeight );
-  QgsDebugMsg( "Entered" );
+  QgsDebugMsg( QString( "thePoint =  %1 %2" ).arg( thePoint.x() ).arg( thePoint.y() ) );
 
   QMap<int, QVariant> results;
 

src/providers/wms/qgswmsprovider.cpp

@@ -3794,18 +3794,31 @@ QMap<int, QVariant> QgsWmsProvider::identify( const QgsPoint & thePoint, Identif
 
   // We don't know highest resolution, so it is difficult to guess any
   // but that is why theExtent, theWidth, theHeight params are here
+  // Keep resolution in both axis equal! Otherwise silly server (like QGIS mapserver)
+  // fail to calculate coordinate because it is using single resolution average!!!
   if ( theWidth == 0 ) theWidth = 1000; // just some number
-  if ( theHeight == 0 ) theHeight = 1000;
+  if ( theHeight == 0 )
+  {
+    theHeight =  myExtent.height() / ( myExtent.width() / theWidth );
+  }
 
   // Point in BBOX/WIDTH/HEIGHT coordinates
   // No need to fiddle with extent origin not covered by layer extent, I believe
   double xRes = myExtent.width() / theWidth;
   double yRes = myExtent.height() / theHeight;
 
+  QgsDebugMsg( "myExtent = " + myExtent.toString() );
+  QgsDebugMsg( QString( "theWidth = %1 theHeight = %2" ).arg( theWidth ).arg( theHeight ) );
+  QgsDebugMsg( QString( "xRes = %1 yRes = %2" ).arg( xRes ).arg( yRes ) );
+
   QgsPoint point;
   point.setX( floor(( thePoint.x() - myExtent.xMinimum() ) / xRes ) );
   point.setY( floor(( myExtent.yMaximum() - thePoint.y() ) / yRes ) );
 
+  QgsDebugMsg( QString( "point = %1 %2" ).arg( point.x() ).arg( point.y() ) );
+
+  QgsDebugMsg( QString( "recalculated orig point (corner) = %1 %2" ).arg( myExtent.xMinimum() + point.x()*xRes ).arg( myExtent.yMaximum() - point.y()*yRes ) );
+
   // Collect which layers to query on
   //according to the WMS spec for 1.3, the order of x - and y - coordinates is inverted for geographical CRS
   bool changeXY = false;