get raster info coordinates messed up

[warmerdam]

The wms_get_feature_info_raster_feature_count.gml test in msautotest/wxs/tomk_workshop.mapis failing in trunk with the following diff:

diff {expected,result}/wms_get_feature_info_raster_feature_count.gml 
13c13
<                                       <gml:coordinates>-81.562500,42.187500 -81.562500,42.187500</gml:coordinates>
---
>                                       <gml:coordinates>443.956163,515.234375 443.956163,515.234375</gml:coordinates>
16,24c16,24
<                       <x>-81.5625</x>
<                       <y>42.1875</y>
<                       <value_0>10</value_0>
<                       <value_1>10</value_1>
<                       <value_2>51</value_2>
<                       <value_list>10,10,51</value_list>
<                       <red>10</red>
<                       <green>10</green>
<                       <blue>51</blue>
---
>                       <x>139.64844</x>
>                       <y>219.14062</y>
>                       <value_0>42</value_0>
>                       <value_1>66</value_1>
>                       <value_2>14</value_2>
>                       <value_list>42,66,14</value_list>
>                       <red>42</red>
>                       <green>66</green>
>                       <blue>14</blue>
29c29
<                                       <gml:coordinates>-80.859375,42.187500 -80.859375,42.187500</gml:coordinates>
---
>                                       <gml:coordinates>445.312500,515.234375 445.312500,515.234375</gml:coordinates>
32,40c32,40
<                       <x>-80.859375</x>
<                       <y>42.1875</y>
<                       <value_0>7</value_0>
<                       <value_1>7</value_1>
<                       <value_2>52</value_2>
<                       <value_list>7,7,52</value_list>
<                       <red>7</red>
<                       <green>7</green>
<                       <blue>52</blue>
---
>                       <x>140.625</x>
>                       <y>219.14062</y>
>                       <value_0>45</value_0>
>                       <value_1>69</value_1>
>                       <value_2>14</value_2>
>                       <value_list>45,69,14</value_list>
>                       <red>45</red>
>                       <green>69</green>
>                       <blue>14</blue>

The gist of the problem seems to be with this calculation:

            /* If projections differ, convert this back into the map  */
            /* projection for distance testing, and comprison to the  */
            /* search shape.  */
            if( needReproject )
                msProjectPoint( &(layer->projection), &(map->projection), 
                                &sPixelLocation );

In cases where the map has non square pixels, and thus "map coordinates" are not real georef coordinates. I'm not sure when this broke.

[warmerdam]

I suspect the problem was introduced in r9716 which added a call to msMapSetFakedExtent() for GetFeatureInfo and relates to #3241.

[warmerdam]

I have emailed Assefa who is responsible for r9716 to explain the requirement. Things are fine for me if I regress and I can't understand the justification for the commit.

[warmerdam]

Assefa has indicated he needed to add the r9716 code in order to get computation of a shplabel value to work properly. I have tried, but not been able to reproduce the problem myself. I have attached a somewhat cutdown version of the msautotest/wxs/tomk_workshop test which runs GetFeatureInfo (see test.sh) producing a result using an html template. Without r9716 (see mapwms_old.c) I get reasonable looking results including:

                <table class="centered" border="1">
                        <tr>
                                <th>Bounding box</th>
                                <td>-179.640000,-179.640000,179.640000,179.640000</td>
                        </tr>
                        <tr>
                                <th>Image Coordinates</th>
                                <td>141, 91</td>
                        </tr>
                        <tr>
                                <th>Map Coordinates</th>
                                <td>-78.120000, 35.100000</td>
                        </tr>
                        <tr>
                                <th>Number of results</th>
                                <td>1</td>
                        </tr>
                        <tr>
                                <th>Scale</th>
                                <td>9478065.979981</td>
                        </tr>
                        <tr>
                                <th>Cellsize</th>
                                <td>1.805427</td>
                        </tr>
                </table>

                <hr/>
		<table class="centered" border="0">
			<tr>
				<td>Query Results for Layer rivers</td>
			</tr>
		</table>
		<hr/>
		<table class="centered" border="1">
			<tr>
				<th>Result</th>
				<th>Name</th>
				<th>System</th>
			</tr>
			<tr>
				<td>1</td>
				<td>Ohio</td>
				<td>Mississippi</td>
				<td>-79.438794,41.524824</td>
			</tr>
                </table>

The last entry above after Mississippi is a [shplabel] substitution. But with the changeset (see mapwms_new.c) I get:

                        <tr>
                                <th>Bounding box</th>
                                <td>-0.000000,-100.000000,500.000000,400.000000</td>
                        </tr>
                        <tr>
                                <th>Image Coordinates</th>
                                <td>141, 91</td>
                        </tr>
                        <tr>
                                <th>Map Coordinates</th>
                                <td>141.282565, 208.695652</td>
                        </tr>
                        <tr>
                                <th>Number of results</th>
                                <td>1</td>
                        </tr>
                        <tr>
                                <th>Scale</th>
                                <td>13190361.250251</td>
                        </tr>
                        <tr>
                                <th>Cellsize</th>
                                <td>2.512563</td>
                        </tr>
                </table>

                <hr/>
		<table class="centered" border="0">
			<tr>
				<td>Query Results for Layer rivers</td>
			</tr>
		</table>
		<hr/>
		<table class="centered" border="1">
			<tr>
				<th>Result</th>
				<th>Name</th>
				<th>System</th>
			</tr>
			<tr>
				<td>1</td>
				<td>Ohio</td>
				<td>Mississippi</td>
				<td>$x,$y</td>
			</tr>

Note that the [shplabel] substitution fails and leaves $x,$y and bounding box and map coordinate values are crazy. Also the cellsize and scale are different though it is hard to determine which is right. One thing I did notice is that without the changeset the query map produced in the tmp directory is wrong, and with the changeset it is correct.

Generally speaking the faked extents are only needed when drawing. So it sort of makes sense that the query map only works with faked extents in effect. But all the stuff that depends on using the map coordinates as projected coordinates break down with fakedextents since in this case the "map coordinate" system isn't really lat/long, but instead just pixel/line coordinates on the map display.

What I'd appreciate is a demonstration of a case where the [shplabel] substitution does not work properly without the changeset applied.

I'd also like to take this opportunity to appologise for the complexity and confusingness of the faked extent approach to handling non-square pixels and rotated maps. It's a mess!

[warmerdam]

demonstration map showing results with/without the changeset.

[assefa]

I am adding the note I sent to Frank here to keep trace:

"I have tried to understand the requirement for this code yesterday night. What I am trying is to return the label position (using shplabel template tag) for queried features I believe that I was having issues when trying to get the position of a label in cases where the image is stretched.

I have attached a zip file that shows what I want to do and the problem I was having.

Basically in this case, there are 2 WMS request (non-square pixels case) done, one is a GetMap and the other is a GetFeatureInfo. The idea is that using the templating in Mapserver with the shplabel tag, It should be possible on the client side to get a suitable x,y position that can be used for labeling. This position should fit nicely on top of the image obtained by thg GetMap request (and should basically correspond to the same label position as if the label was drawn directly on the server side) The way the position is calculated is (maptemplate.c:1559):

• we convert the queried shape from layer to map projection
• we clip the shape to the map extents
• convert in the shape in pixel coordinates
• compute the text position
• reconvert back to the position to the output projection

This way of calculating is what is done when MapServer renders the text.

Without the msMapSetFakedExtent, I was not getting the correct position I was expecting.

In the examples I sent, doc.kml is with the code as It stands now and doc-problem.kml is without the problematic code. If you load both in Google earth, and zoom in a bit, you will see line features and a text coming from a WMS GetMap request and on top of that you will see the pin points (lable positions) returned by the GetFeatureInfo request (through the shplabel template). In one case the pin position is exactly on top of the text returned by the GetMap and in the other case, It is slightly off.

I hope these examples show what I am trying to do.

I am not sure I took the correct approach in adding that part of the code at that level but before removing it, I would like to have your feed backs on it. "

[assefa]

kml demo explaining the problem

[assefa]

I will try to get a way to demonstrate this. But from my understanding the case I was trying to solve with the problematic piece of code is that if I want to superpose a an image obtained by a GetMap request (non square pixels which "stretches" the image), and a point position obtained by using a shplabel tag through a GetFeature request (the same exact wms parameters as the GetMap) I was not getting the position I was expecting. I believe I would not have run into this issue if the WMS request was not "stretching" the image.

[warmerdam]

This issue is also messing up a WMS GetFeatureInfo request in the msautotest/wxs/wfs_ogr.map test set - points are being returned in the wrong projection. I'll spend some time tonight to review what Assefa has written and how the code works, but I don't think I can leave the r9716 changes in as they are now since they break several things.

[warmerdam]

I looked into it a bit tonight, but the amazon instances are not responding so I can't try the kml's (and I don't normally use google earth).

Perhaps I'll look again tomorrow.

[assefa]

Do you see a problem if I put the code inside #ifdefs? The ec2 instances are down and I can't have them up until at least Monday. Once I have a clear demo of the problem I was trying to solve, I can reengage you through this bug. I feel that I am blocking things and not being efficient until I really allocate a good chunk of time to work on things.

[warmerdam]

Having it in ifdef's is fine, and I can certainly wait a few days. But I do feel strongly that this needs to be resolved before the MapServer release which I think we are trying to move towards reasonably soon.

[assefa]

I have tried to reproduce the problem I was trying to solve with this piece of code for the last few hours but I was unable to do so. With or without the code I get consistently the same result. I will give it again a try later tonight and report. If I can not justify it after that I will remove the code.

[assefa]

I have removed the problematic code: I am getting the correct results for polygons and points when superposing "stretched" wms images produced with GetMap and a feature returned through a GetFeature request (and using templating). I am still a bit off in Google earth for the line features (regardless if this code is present or not). I am using mapserver-svn. Committed in r 10700