segfault while processing GetFeatureInfo

[tbonfort]

I have a raster layer (actually, a bunch of them) where MapServer sometimes segfaults when serving a GetFeatureInfo request, depending on the coordinate of the point being queried. I have come up with a couple of test cases, which are exactly identical to each other except for the value of X in the GetFeatureInfo request. In one case, X=913, and MapServer correctly processes the GetFeatureInfo request and returns correct results. In the other case, X=912, and MapServer segfaults. Both of these values of X are valid, and both correspond to pixels in the image where MapServer correctly displays the raster in question when called with a GetMap request.

My hunch is that the problem might be related to a misconfiguration of MapServer on my part, or a misconfiguration of the projections involved. The layer data is in a LAEA projection, and due to constraints on the project that are beyond my control, I have MapServer reprojecting it to a bizarre web mercator projection (EPSG 3857, aka EPSG 3785) for output, so the GetFeatureInfo request is expressed in EPSG 3857. But I'm stumped as to what the problem could be, and having MapServer crash with a segmentation fault seems like an indication of some kind of bug, even if I've got a configuration error.

I've tried turning on debugging in the mapfile, but nothing useful shows up in the log file.

I've tested this with several versions of MapServer, including 6.0.3 and 6.2-beta1, with the same results. I'm running these tests on a CentOS 6.2 system, compiling MapServer with the options --with-proj --with-ogr --with-gdal --with-wfs --with-wcs --with-wms --with-wmsclient --with-wfsclient and --with-php=/usr/include/php, and using gdal 1.7.2 and
proj 4.7.0.

My mapfile is appended below, along with the two test requests I mentioned above.

The layer data file that I am using is too big to attach to this email, but I've created a tgz file containing everything needed to run my tests -- the mapfile, layer file, and test script; you can download it here:

  http://rain.nemac.org/~mbp/mapserver-gfi-problem.tgz

Here is the mapfile I am using:

MAP
  CONFIG "MS_ERRORFILE" "mapserver.log"
  DEBUG 5
  PROJECTION
    "init=epsg:3857"
  END
  WEB
    HEADER "./dummy_template"
    FOOTER "./dummy_template"
    METADATA
      "ows_enable_request"     "*"
      "wms_srs"                "EPSG:3857"
    END
  END
  LAYER
    NAME layer1
    PROJECTION
      "+proj=laea"
      "+lat_0=45"
      "+lon_0=-100"
      "+x_0=0"
      "+y_0=0"
      "+a=6370997"
      "+b=6370997"
      "+units=m"
      "+no_defs"
    END
    TYPE     RASTER
    DUMP     TRUE
    STATUS   OFF
    DATA     layer1.tif
    HEADER   ./dummy_template
    TEMPLATE ./dummy_template
    METADATA
      "wms_title"             "layer1"
      "wms_abstract"          "layer1"
      "gml_include_items"     "value_0"
    END
  END
END

And here are the two query strings for the two test cases:

this one runs correctly:
map=./mapfile.map&TRANSPARENT=true&SERVICE=WMS&VERSION=1.3.0&REQUEST=GetFeatureInfo&STYLES=&FORMAT=image/png&INFO_FORMAT=application/vnd.ogc.gml&BBOX=-11859604.53136514,4696887.582602486,-11830673.11615924,4710053.860724591&CRS=EPSG:3857&LAYERS=layer1&QUERY_LAYERS=layer1&WIDTH=1514&HEIGHT=689&X=913&Y=334

this one causes a seg fault:
map=./mapfile.map&TRANSPARENT=true&SERVICE=WMS&VERSION=1.3.0&REQUEST=GetFeatureInfo&STYLES=&FORMAT=image/png&INFO_FORMAT=application/vnd.ogc.gml&BBOX=-11859604.53136514,4696887.582602486,-11830673.11615924,4710053.860724591&CRS=EPSG:3857&LAYERS=layer1&QUERY_LAYERS=layer1&WIDTH=1514&HEIGHT=689&X=912&Y=334

[tbonfort]

hi,
first of all, congrats for the perfect test case!

I've tracked down your issue and managed to reproduce it. For some reason or another, proj.4 is failing to reproject the query point, and the featureinfo writing code chokes on this, which is clearly a bug. We can fix the segfault but the resultset will be empty which is incorrect.

I'm not a gis/projection guru, can you confirm that the queried point is actually valid for your raster projection? If so there's a bug in proj that would need to be adressed.

[tbonfort]

I've just proposed a fix that also adds a more verbose proj error message, which gives:

Starting program: /home/tbonfort/local/bin/mapserv QUERY_STRING='map=./mapfile.map&TRANSPARENT=true&SERVICE=WMS&VERSION=1.3.0&REQUEST=GetFeatureInfo&STYLES=&FORMAT=image/png&INFO_FORMAT=application/vnd.ogc.gml&BBOX=-11859604.53136514,4696887.582602486,-11830673.11615924,4710053.860724591&CRS=EPSG:3857&LAYERS=layer1&QUERY_LAYERS=layer1&WIDTH=1514&HEIGHT=689&X=912&Y=334'
[Thread debugging using libthread_db enabled]
Using host libthread_db library "/lib/x86_64-linux-gnu/libthread_db.so.1".
[Mon Jul 23 12:13:25 2012].396349 msRasterQueryByRect(layer1): entering.
[Mon Jul 23 12:13:25 2012].401118 msProjectRect(): some points failed to reproject, doing internal sampling.
Content-type: application/vnd.ogc.gml

<?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">
    <layer1_layer>
    <gml:name>layer1</gml:name>
[Mon Jul 23 12:13:25 2012].402457 msProjectPoint(): Projection library error. proj says: tolerance condition error
<!-- Warning: Failed to reproject point: msProjectPoint(): Projection library error. proj says: tolerance condition error -->
    </layer1_layer>
</msGMLOutput>
[Mon Jul 23 12:13:25 2012].402554 mapserv request processing time (msLoadMap not incl.): 0.006s
[Mon Jul 23 12:13:25 2012].402573 msFreeMap(): freeing map at 0x6051b0.
[Mon Jul 23 12:13:25 2012].402599 freeLayer(): freeing layer at 0x60d2e0.

[embeepea]

Hi Thomas,

Thanks for looking into this!

I am pretty sure that the query point is valid, because I generated it by
clicking on an image of this layer that was generated by a GetMap request
using the same mapfile. I will check with my local GIS experts, though, to
make sure, and will let you know.

In the meantime, it would be helpful to me to know the exact coordinates,
in the query CRS, of the point that mapserver is giving to proj4 for
conversion. I mean, I'm assuming that mapserver is converting the passed X
and Y (which are pixel coordinates, relative to the passed WIDTH and HEIGHT
values) to the coordinates of a point in EPSG:3857 (the query coordinate
reference system), and then handing THOSE coordinates off to proj4 for
conversion to the underlying layer coordinate reference system. So proj4
is being asked to convert a point in EPSG:3857 to the LAEA projection that
the layer data is in.

Doing the math myself, I'm thinking that the EPSG:3857 coordinates of the
query point would be (-11842176.888916,4703270.074464). But could you
confirm that for me with your copy of mapserver (temporarily put in a print
statement to print them out, or examine them with a debugger)?

Also, it might be a good idea to add those coordinates to the error message
that mapserver generates in your proposed fix; it could be helpful to
people trying to debug a problem.

Thanks again,

--Mark

On Mon, Jul 23, 2012 at 6:00 AM, Thomas Bonfort <
reply@reply.github.com

wrote:

hi,
first of all, congrats for the perfect test case!

I've tracked down your issue and managed to reproduce it. For some reason
or another, proj.4 is failing to reproject the query point, and the
featureinfo writing code chokes on this, which is clearly a bug. We can fix
the segfault but the resultset will be empty which is incorrect.

I'm not a gis/projection guru, can you confirm that the queried point is
actually valid for your raster projection? If so there's a bug in proj that
would need to be adressed.

---

Reply to this email directly or view it on GitHub:
#4403 (comment)

[tbonfort]

I suspect there's some adjustment related to the centering of the query point w.r.t. the pixel. The query point that's failing is
-11842309.595810996,4703658.4663789356

[embeepea]

I agree; those coords are close enough to what I got that I imagine the
difference is in the pixel adjustment.

Note, though, that the proj command-line utility seems perfectly happy to
convert that point, and as far as I can tell, it's getting a reasonable
answer ---- at least it's not crashing. Here's a shell script that does
the conversion:

#! /bin/bash

echo '-11842309.595810996 4703658.4663789356' | invproj -f "%6f" +proj=merc
+a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0
+units=m +nadgrids=@null +wktext +no_defs

echo '-106.381277 38.874135' | proj +proj=laea +lat_0=45 +lon_0=-100 +x_0=0
+y_0=0 +a=6370997 +b=6370997 +units=m +no_defs

The first call, to invproj, converts the original point from EPSG:3857 to
lon/lat; the arguments to invproj for specifying the projection are from
the line starting with <3857> in my /usr/local/share/proj/epsg.

The second call, to proj, converts the lon/lat returned by the first call
to the LAEA projection of layer1.tif. The output of the script is:

-106.381277 38.874135
-552535.42 -659643.78

Just as an extra sanity check, I even did the conversion back in the other
direction, to convert the final LAEA coordinates back to EPSG:3857, and
proj correctly gives the original point as the result.

I'm able to do all this with the same version of proj that has been
compiled into my copy of mapserver. The only difference is that I'm using
the proj commmand-line utility, rather than calling the proj library from a
C program.

So, I'm wondering if the problem has something to do with the way that
mapserver is calling the proj library, or with the way it's loading one of
the projections involved.

I'm happy to help debug this, but could you help me out by pointing me to
the line in the mapserver source code where mapserver is calling the proj
library to do the conversion? I took a look a look at the mapserver source
yesterday and it wasn't immediately obvious to me.

Thanks!

--Mark

On Tue, Jul 24, 2012 at 4:04 AM, Thomas Bonfort <
reply@reply.github.com

wrote:

I suspect there's some adjustment related to the centering of the query
point w.r.t. the pixel. The query point that's failing is
-11842309.595810996,4703658.4663789356

---

Reply to this email directly or view it on GitHub:
#4403 (comment)

[tbonfort]

I'll be offline till the beginning of next week so won't be able to help much more until then. here's the callstack of the failing proj call:

s_inverse (xy={x = -1.8587843622922748, y = 0.73829236874211923}, P=0x60e880)
pj_inv (xy={x = -1.8587843622922748, y = 0.73829236874211923}, P=0x60e880)
pj_transform (srcdefn=0x60e880, dstdefn=0x60cc20, point_count=1, point_offset=1, x=0x63b8f0, y=0x63b8f8, z=0x7fffffffdbd8)
msProjectPoint (in=0x60d408, out=0x6068a0, point=0x63b8f0)
msProjectLine (in=0x60d408, out=0x6068a0, line=0x63c8a0)
msProjectShape (in=0x60d408, out=0x6068a0, shape=0x7fffffffe110)
msGMLWriteQuery (map=0x6051b0, filename=0x0, namespaces=0x7ffff7b9f184 \"MGO\")
msWMSFeatureInfo (map=0x6051b0, nVersion=66304, names=0x603d50, values=0x604080, numentries=16, wms_exception_format=0x0, ows_request=0x7fffffffe440)
msWMSDispatch (map=0x6051b0, req=0x603d10, ows_request=0x7fffffffe440, force_wms_mode=0)
msOWSDispatch (map=0x6051b0, request=0x603d10, ows_mode=-1)
msCGIDispatchRequest (mapserv=0x603ad0)
main (argc=2, argv=0x7fffffffe638)

the tolerance error seems to be coming from:

rh = hypot(xy.x, xy.y);
if ((lp.phi = rh * .5 ) > 1.) I_ERROR;

which is around line 148 in PJ_laea.c of a more-or-less trunk version of proj

[tbonfort]

@embeepea any update on this one ?

[embeepea]

Sorry, I haven't had a chance to get back to this --- I've been busy with
another project. I will be out of town until the middle of next week and
probably won't get a chance to look at it again until then, but I will get
back to it then.
--Mark

On Tue, Jul 31, 2012 at 6:10 AM, Thomas Bonfort <
reply@reply.github.com

wrote:

@embeepea any update on this one ?

---

Reply to this email directly or view it on GitHub:
#4403 (comment)

[tbonfort]

I've committed the fix that prevents the segfault into 6.2. Leaving this issue opened but demilestoning it for now.

[embeepea]

OK, I finally got back to this after several weeks on another project.

The short version of the situation is that I think I may have found the problem, and I have a proposed fix for it. This involves code that I only have a very sketchy understanding of, though, so please scrutinize my solution closely to make sure it looks correct, and let me know if I'm way off base!

I've forked the MapServer repo and pushed a commit to my fork that has my changes; you can check it out at

github.com/embeepea/mapserver.git

Both of the tests that I submitted originally with this issue run correctly with this fix, and it seems to correctly fix the coordinate issue I describe below.

Here's a longer explanation:

It looks to me like the problem has to do with MapServer's reprojection of the coordinates of the "hit" pixel into the SRS of the map. In particular, look at the following results. For the query string

map=./mapfile.map&TRANSPARENT=true&SERVICE=WMS&VERSION=1.3.0&REQUEST=GetFeatureInfo&STYLES=&FORMAT=image/png&INFO_FORMAT=application/vnd.ogc.gml&BBOX=-11859604.53136514,4696887.582602486,-11830673.11615924,4710053.860724591&CRS=EPSG:3857&LAYERS=layer1&QUERY_LAYERS=layer1&WIDTH=1514&HEIGHT=689&X=913&Y=334

MapServer returns the following GML:

<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">
    <layer1_layer>
    <gml:name>layer1</gml:name>
        <layer1_feature>
            <gml:boundedBy>
                <gml:Box srsName="EPSG:3857">
                    <gml:coordinates>8935763.864626,-5609465.555845 8935763.864626,-5609465.555845</gml:coordinates>
                </gml:Box>
            </gml:boundedBy>
            <value_0>126</value_0>
        </layer1_feature>
    </layer1_layer>
</msGMLOutput>

This is the first test case of the two cases I submitted --- this is the one that does NOT crash.

Look closely at coordinates in the query string. If I'm interpreting it correctly, the query is asking MapServer for information about what is at the location X=913 Y=334 in a 1514 x 689 pixel image whose extent (BBOX) is

-11859604.53136514   < x <  -11830673.11615924
  4696887.582602486  < y <   4710053.860724591

expressed in the SRS EPSG:3857. The response says that it found a result in layer1, whose pixel value is 126, and whose coordinates in SRS EPSG:3857 are X=8935763.864626, Y=-5609465.555845. But those coordinates are not even in the BBOX of the query!

I'll admit that I'm making a few unverified (I haven't read the WMS GetFeatureInfo spec) assumptions here about the meaning of the values in the WMS request and response, and if I'm wrong about this and need to formulate my query string differently, just let me know, and I apologize for taking up your time! But if my interpretation is correct, I think there's a problem with MapServer's handling of the coordinates of the "hit" region, and this problem happens universally --- not just in the version of my test that was causing a segfault. I suspect that the reason it segfaults only in some cases and not others is that the values being handed to proj cross certain thresholds in some cases but not others. But in all cases the result that proj returns is invalid.

After studying the MapServer source for quite a while, it looks to me like what is happening is that sRasterQueryByRectLow() in maprasterquery.c is converting the query point into the mp SRS in order to do its distance calculation in those coordinates. If it finds a hit, it loads the query results with those query point coordinates which have been converted to the map SRS. The problem is that the code that formulates the GML response expects the results to be in the layer's SRS, and it's invoking proj to convert them to the map SRS. So the problem seems to be one too many coordinate transformations!

The fix that I came up with is to change msRasterQueryByRectLow() so that it loads the results cache with the original layer SRS version of the query point.

I'm not familiar enough with all the assumptions made in the code to know for sure that this is the correct solution, but it seems to work for my test case. Please take a look at it and let me know what you think!

Thanks again,

--Mark

[tbonfort]

@embeepea I'm not enough into raster queries to be able to make an informed comment on your fix. maybe @warmerdam can have a look ?

[sdlime]

My apologies, when I saw Thomas fully engaged I let this drop of my radar. If Frank is around he's obviously in the best place to comment. Otherwise I know the query operations well. With the test case and the repo I should be able to work on this one. --Steve

[embeepea]

Thanks. Steve or Frank, see my comment #4403 (comment) from a couple of weeks ago for the details of what I found when I looked into this.

A short summary of the situation is this: it looks to me, at least in the test case that I came up with, that MapServer is converting the coordinates of the raster query point from the layer SRS to the map SRS twice:

• (a) once inside the raster query code (in function msRasterQueryByRectLow), and
• (b) again in the code that formats the GetFeatureInfo response into GML (function msGMLWriteQuery)

This results in invalid coordinates, and sometimes causes proj to barf.

The fix that I came up with (see my forked repo at https://github.com/embeepea/mapserver) was to modify msRasterQueryByRectLow so that it returns the hit point coordinates in the layer SRS, instead of converting to the map SRS. This seems to fix my test case, but I'm not familiar enough with the assumptions in the MapServer code to know for sure that it's the correct thing to do. The other possibility would be to let msRasterQueryByRectLow() do the conversion, and change msGMLWriteQuery to not do it; which one is correct just depends on what assumptions are in place regarding how data should be returned from the low level query code.

If I can be of more help in figuring this out, let me know!

Thanks,

--Mark

On Fri, Sep 7, 2012 at 9:57 AM, Steve Lime notifications@github.com wrote:

My apologies, when I saw Thomas fully engaged I let this drop of my radar.
If Frank is around he's obviously in the best place to comment. Otherwise I
know the query operations well. With the test case and the repo I should be
able to work on this one. --Steve

—
Reply to this email directly or view it on GitHubhttps://github.com//issues/4403#issuecomment-8365919.

[tbonfort]

For the record, the proposed modification is in aa59fd5

[tbonfort]

@sdlime, seems like Frank won't be looking into this. The proposed change seems reasonable to me, could this make it into 6.2 ? I'll take care of running the testsuite to see if this has some other side-effects, but would be grateful if you could have a quick look. thanks, thomas

[tbonfort]

ok, I'll bite the bullet and bring this in for 6.2. I've added two raster querying tests that clearly show that something was wrong. In https://github.com/mapserver/msautotest/blob/86960175ace0601afe9f57630200066afce50139/wxs/expected/wms_rast_featureinfo_reproj.xml#L14 , before this patch the bounds would end up as

<gml:coordinates>-121.48979,0.00040 -121.48979,0.00040</gml:coordinates>

instead of now as

<gml:coordinates>-117.44501,44.71061 -117.44501,44.71061</gml:coordinates>

I still needed to fix one last thing: in https://github.com/mapserver/msautotest/blob/86960175ace0601afe9f57630200066afce50139/wxs/expected/wms_rast_featureinfo_reproj.xml#L17 , the x and y values of the query point where being returned in map srs, they are now in layer srs. In my previous commit, I have made these modifications:

• renamed variables to make it clearer that we are working on reprojected pixel coordinates
• added entries to the raster resultsetinfo to contain the reprojected pixel values, so that those can be added to the shape x,y values instead of the ones in layer SRS.

[sdlime]

I can have a look later today if that helps... --Steve

[tbonfort]

@sdlime just a quick look should be enough. The raster querying code is a bit particular in the sense that it puts the geometry information in the result shape and values. One goes through reprojection, and of course not the other, which ends up in us having to store both the projected and original pixel locations.

[warmerdam]

Guys,

Sorry - I didn't notice this was associated with me, and I am not clear on
how to interact with the new bug system. Thanks for addressing this.

Best regards,
Frank
On Sep 29, 2012 7:15 AM, "Thomas Bonfort" notifications@github.com wrote:

@sdlime https://github.com/sdlime just a quick look should be enough.
The raster querying code is a bit particular in the sense that it puts the
geometry information in the result shape and values. One goes through
reprojection, and of course not the other, which ends up in us having to
store both the projected and original pixel locations.

—
Reply to this email directly or view it on GitHubhttps://github.com//pull/4403#issuecomment-9004368.

[tbonfort]

@warmerdam, thanks. The proposed patch can be seen here: https://github.com/tbonfort/mapserver/commit/c65723a45d2c979072f002a58c054a8c5ca6b162 , it has the inconvenience of having to add two members to the lrinfo struct, but I don't see how differently we can be doing this without being backwards incompatible with the current behavior that adds x and y to the feature attributes.
regards, thomas

[embeepea]

Guys,

Thanks for looking into this; let me know if I can be of any more help with
it.

--Mark

On Sat, Sep 29, 2012 at 11:49 AM, Thomas Bonfort
notifications@github.comwrote:

@warmerdam https://github.com/warmerdam, thanks. The proposed patch can
be seen here: tbonfort/mapserver@c65723ahttps://github.com/tbonfort/mapserver/commit/c65723a45d2c979072f002a58c054a8c5ca6b162, it has the inconvenience of having to add two members to the lrinfo
struct, but I don't see how differently we can be doing this without being
backwards incompatible with the current behavior that adds x and y to the
feature attributes.
regards, thomas

—
Reply to this email directly or view it on GitHubhttps://github.com//pull/4403#issuecomment-9005200.