Raster Query

Raster query support is/was incorporated into MapServer in April/2004 (post MapServer 4.2) by Frank Warmerdam.

The official reference documentation on raster query is in the Raster Data Access document. This Trac page is outdated but preserved for possibly useful additional details.

General Approach

Raster query support is setup to return a set of request result "shapes" for each pixel matching the query. In MapServer CGI, normal template substitution is done. In MapScript the normal mechanisms for reading a query result can be applied.

Currently QueryByPoint (with single or multiple modes), QueryByRect and QueryByShape are supported. Care should be taken to avoid providing a large query area (selecting alot of pixels) as each selected pixel requires over 100 bytes of memory for temporary caching. The RASTER_QUERY_MAX_RESULT PROCESSING item can be used to restrict the maximum number of query results that will be returned. The default is one million which would take on the order of 100MB of RAM.

Each pixel matching the query has (potentially) the following attributes:

  • x: X location of pixel
  • y: Y location of pixel
  • value_list: a comma seperated list of the values of all selected bands at the pixel (raw values).
  • value_[n]: the value for the n'th band in the selected list at this pixel (zero based).
  • class: Name of the class this pixel is a member of. (zero based class number for unnamed classes).
  • red: red component of display color for this pixel.
  • green: green component of display color for this pixel.
  • blue: blue component of display color for this pixel.

MapScript Example

This Python MapScript example gives an idea of how the query can be done, and results fetched in MapScript.

map = mapscript.Map('rquery.map')
layer = map.getLayer(0)

pnt = mapscript.Point()
pnt.x = 440780
pnt.y = 3751260
    
layer.queryByPoint( map, pnt, mapscript.MS_MULTIPLE, 180.0 )

layer.open()
for i in range(1000):
    result = layer.getResult( i )
    if result is None:
        break
        
    print '(%d,%d)' % (result.shapeindex, result.tileindex)
        
    s = layer.getShape( result.shapeindex, result.tileindex )
    for i in range(layer.numitems):
        print '%s: %s' % (layer.getItem(i), s.getValue(i))
            
layer.close()

Note that the normal vector feature layer mechanism (layer.open(), layer.getShape()) is used to fetch the query results. The query operation creates an internal fake set of shapes representing the pixels that matched the query. Attempts to use layer.open() and layer.getShape() on a raster layer that has not been previously queried will just return zero shapes (and an empty item list). Furthermore, saving query results and re-loading them won't work since the way MapServer saves query results is to just save the id's, not the actual contents.

Query With Templates

With query templates pixel "attributes" can be substituted normally. For instance the following template (specified in the LAYER with the TEMPLATE keyword) would work normally:

 Pixel:<br>
  values=[value_list]<br>
  value_0=[value_0]<br>
  value_1=[value_1]<br>
  value_2=[value_2]<br>
  RGB = [red],[green],[blue]<p>
  Class = [class]<br>

Working Notes

  • Note that WMS GetFeatureInfo requests work against raster layers!
  • Currently scaling is not applied to non-eight bit values in before computing class or color. Thus class, red, green and blue results may be a little or alot wrong in some circumstances. See
  • The _msQueryByIndex(), msQueryByIndexAdd(), msQueryByAttributes() and msQueryByFeatures() are not supported for rasters. I don't really understand what they are or what they do.
  • The raster layer follows similar rules to other layers to enable query. Basically a template must be enabled on the layer. Also, min/max scale, and status checks are applied.
  • There is a bug with classified layers that use [red], [green] and [blue] in expression. See bug 1021 ( http://trac.osgeo.org/mapserver/ticket/1021).