Query (GetFeatureInfo) gives incorrect results (same result, bbox nan)

[bartvde@…]

Using 4.10.

When using Mapserver native SHAPEFILE connector, we get incorrect results on
doing a WMS GetFeatureInfo request. If we using CONNECTIONTYPE OGR all works
fine. It always gives back the same feature, and the boundingbox has nan values.

Example request:
http://145.50.148.45:8081/cgi-bin/mapserv410?map=/data/OGC_UMN_services/nis/getfeatureinfobug.map&LAYERS=test&REQUEST=GetFeatureInfo&SERVICE=WMS&FORMAT=image%2Fpng&STYLES=default&HEIGHT=450&QUERY_LAYERS=test&VERSION=1.1.1&SRS=EPSG%3A28992&WIDTH=750&SLD=http%3A%2F%2Ftest.intranet.rijkswaterstaat.nl%2Fapps%2Fgeoservices%2Fnis%2Fsld%2Fnis-mjpv.sld&BBOX=-131666.662034784%2C300000.00277913%2C451666.662034784%2C649999.99722087&TRANSPARENT=TRUE&INFO_FORMAT=application%2Fvnd.ogc.gml&X=492&EXCEPTIONS=application%2Fvnd.ogc.se_xml&Y=125

Response using Mapserver SHAPEFILE connector:
<?xml version="1.0" encoding="ISO-8859-1"?>

<msGMLOutput 
	 xmlns:gml="http://www.opengis.net/gml"
	 xmlns:xlink="http://www.w3.org/1999/xlink"
	 xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
	<test_layer>
		<test_feature>
			<gml:boundedBy>
				<gml:Box srsName="EPSG:28992">
					<gml:coordinates>nan,nan nan,nan</gml:coordinates>
				</gml:Box>
			</gml:boundedBy>
			<DIENSTCODE>RWS ZL</DIENSTCODE>
			<DISTRICTCO>510</DISTRICTCO>
			<WEGNUMMER>58</WEGNUMMER>
		</test_feature>
	</test_layer>
</msGMLOutput>



Response using CONNECTIONTYPE OGR:
<?xml version="1.0" encoding="ISO-8859-1"?>

<msGMLOutput 
	 xmlns:gml="http://www.opengis.net/gml"
	 xmlns:xlink="http://www.w3.org/1999/xlink"
	 xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
	<test_layer>
		<test_feature>
			<gml:boundedBy>
				<gml:Box srsName="EPSG:28992">
					<gml:coordinates>248848.347085,552271.298354
248866.128423,552369.886861</gml:coordinates>
				</gml:Box>
			</gml:boundedBy>
			<DIENSTCODE>RWS NN</DIENSTCODE>
			<DISTRICTCO>716</DISTRICTCO>
			<WEGNUMMER>34</WEGNUMMER>
		</test_feature>
	</test_layer>
</msGMLOutput>

Will attach MAP file and shapefile to reproduce.

[bartvde@…]

MAP file plus shapefile in tgz format

[sdlime]

CC'ing

[sdlime]

Did this work in previous versions?

Steve

[bartvde@…]

Steve, the same problem occurs at least in 4.8. I did not check any other versions.

[sdlime]

Ok, that helps. At least we know we didn't bust anything and it was already 
broken... I'm wondering if it might be related to the plugin architecture work, 
a missing or broken getExtent function.

Steve

[bartvde@…]

It's not only the bbox Steve, a totally wrong feature is returned not
corresponding to the X and Y values, and it is always the same feature
irrespective of the X and Y given.

[sdlime]

Working on it. I can confirm the behavior here. I can also confirm that the
query doesn't seem to work outside of WMS using a URL like:

http://localhost/cgi-bin/mapserv?map=/Users/sdlime/dev/tmp/1931/getfeatureinfobug.map&mapheight=450&qlayer=test&MAPWIDTH=750&imgext=-131666.662034784+300000.00277913+451666.662034784+649999.99722087&img.x=492&img.y=125&mode=query

Leads me to believe there might be a problem or inconsistency with the shapefile
that OGR is more robust in dealing with.

Steve

[sdlime]

Ok, I wrote a little mapscript application to loop though the test shapefile and
output the bounds for each shape. A significant number of them have corrupt
boundary values (according to the native shapefile driver in MapServer). That
driver uses bbox for all query pre-filtering so that's likely the problem. There
may be an error with the native driver code so I've cc'd Frank and will post the
script output.

Steve

[sdlime]

MapScript output

[bartvde@…]

Frank, are you able to have a look at this? 

This is a pretty annoying bug with Mapserver native shapefile support that we
have on many of our shapefiles.

OGR does this right.

Thanks in advance.

[fwarmerdam]

Based on a quick scan of the file with shpdump it seems there are 35
features with no vertices.  These all have nan for bounds. 

eg.

Shape:1375 (ArcM)  nVertices=0, nParts=0
  Bounds:(         nan,         nan, 0, -1.79769e+308)
      to (         nan,         nan, 0, -1.79769e+308)

I don't see any features with actual vertices and corrupt bounds. 

It seems quite plausible that ArcGIS considers nan a good way to mark
shapes as having no meaningful bounds.  

If this seems like a serious problem, you could alter the MapServer
shapefile code to reset the bounds to all zero if a shape has
no parts and no vertices.  

Note that OGR does not keep the bounds from the shapefile with the 
geomtery.  When things are turned into a shapeObj I suppose the bounds
must be recomputed from the geometry.

[bartvde@…]

Steve, Frank,

we get a lot of shapefiles with these characteristics, and they all come from
ESRI environments.

Steve, can anything be done about this in the Mapserver query code? For instance
like Frank suggests resetting the bounds to 0?

Thanks in advance.

[sdlime]

We could augment the shapefile reader to test for NaN (how does one do that in 
C?) and set the bounds according as Frank suggests. However, we need to be able 
to deal with empty shapes anyway, perhaps an isShapeValid macro that checks the 
number of parts. These shapefiles draw ok?

Steve

[bartvde@…]

Hi Steve,

the shapefiles draw okay. I can't answer the other C questions though :-)

Shapelib seems to have something for it though if I interpret the code correctly
(shapelib-1.2.10\contrib\my_nan.h).

Bart

[fwarmerdam]

Steve,

The following is what I use for a NAN test in GDAL.  It should
be reasonably portable.  But I'd advise just clearing the bounds
if the parts==vertices==0 as that seems to cover off the problem
features.  

/* -------------------------------------------------------------------- */
/*      Handle isnan() and isinf().  Note that isinf() and isnan()      */
/*      are supposed to be macros according to C99.  Some systems       */
/*      (ie. Tru64) don't have isinf() at all, so if the macro is       */
/*      not defined we just assume nothing is infinite.  This may       */
/*      mean we have no real CPLIsInf() on systems with an isinf()      */
/*      function but no corresponding macro, but I can live with        */
/*      that since it isn't that important a test.                      */
/* -------------------------------------------------------------------- */
#ifdef _MSC_VER
#  define CPLIsNan(x) _isnan(x)
#  define CPLIsInf(x) (!_isnan(x) && !_finite(x))
#  define CPLIsFinite(x) _finite(x)
#else
#  define CPLIsNan(x) isnan(x)
#  ifdef isinf 
#    define CPLIsInf(x) isinf(x)
#    define CPLIsFinite(x) (!isnan(x) && !isinf(x))
#  else
#    define CPLIsInf(x)    FALSE
#    define CPLIsFinite(x) (!isnan(x))
#  endif
#endif

[bartvde@…]

Steve,

can this be done for 5.0?

Bart

[sdlime]

Yes, will try to get into CVS today. One question for Frank. You mentioned:

"But I'd advise just clearing the bounds if the parts==vertices==0 as that seems
to cover off the problem features."

By clearing the bounds you mean setting to 0.0, correct? 

Steve

[sdlime]

Ok, actually looking at code. The initial bounding box stuff is done without 
looking at the shape itself (for speed) in msSHPReadBounds so I think we'll 
actually need to do the NaN test in that function and set bounds to all 0.0. 
(We do that for NULL shapes as it is.)

We also read the bbox in msSHPReadShape, but in that case we do have the number 
of parts and/or vertices and can test as Frank suggests.

Bart, can you test 5.0?

Steve

[bartvde@…]

Hi Steve, yes we can do a test at our workplace with CVS head (aka 5.0). Let me
know when the change has been committed (and ideally in which file). Thanks in
advance!

[sdlime]

Ok, I believe this is fixed in CVS head. Bart's data works just fine now with
point queries. I only updated the msSHPReadBounds function to check for NaN
bounds and didn't mess with msSHPReadShape or msSHPReadPoint since they seem to
be working just fine. I also took the liberty of cleaning up the return values
for that function to MS_SUCCESS/MS_FAILURE instead of the odd 0/-1 (updated any
calls to the function elsewhere accordingly).

Marking as fixed. Bart, can you confirm please?

Steve