[GIS] How does mapscript queryByShape work

mapscriptmapserveropenlayers-2

I am trying to return the results of queryByShape call from a (large) raster file given a user specified polygon.

enter image description here

For example, when a user plots the following polygon on an openlayers layer as shown in this diagram, I send the coordinates of the polygon vertices to the python mapscript backend.

With python mapscript, I re-create the polygon as a mapscript shapeObj.

line = mapscript.lineObj()
line.add( mapscript.pointObj( 16.484375, 59.736328125 ) )
line.add( mapscript.pointObj( 17.4951171875, 55.1220703125 ) )
line.add( mapscript.pointObj( 24.74609375, 55.0341796875) )
line.add( mapscript.pointObj( 22.5927734375, 61.142578125) )
line.add( mapscript.pointObj( 16.484375, 59.736328125) )
poly = mapscript.shapeObj( mapscript.MS_SHAPE_POLYGON )
poly.add( line )

Then, I experiment with the queryByShape call like this:

mapfile = '/path/to/my/mapfile.map'
mymap = mapscript.mapObj( mapfile )

layer = mymap.getLayer(0)    
layer.queryByShape( mymap, poly )

The returned result is an integer "1", which isn't what I was expecting.

Unfortunately, there isn't much by way of documentation for queryByShape so I am not sure what I am doing wrongly (http://mapserver.org/mapscript/querying.html#by-shape)

How does mapscript queryByShape work? I am actually trying to return a shape object where I can extract information about the selected shape (polygon defined by user) and determine

  1. The total area of this shape,
  2. The % in this shape which is of a specific RGB ("blue") value,
  3. The % in this shape which is of a specific RGB ("green") value etc.

Here's the LAYER definition in my mapfile.

LAYER
    NAME 'GLOBCOVER_L4_200901_200912_V2.3.color'
    TYPE RASTER
    DUMP true
    TEMPLATE 'fooOnlyForWMSGetFeatureInfo'
    EXTENT -189.001389 -115.106069 188.998611 140.108847
    DATA '/media/GLOBCOVER_L4_200901_200912_V2.3.color.tif'

    METADATA
        'ows_title' 'GLOBCOVER_L4_200901_200912_V2.3.color'
    END

    STATUS OFF
    TRANSPARENCY 100

    PROJECTION
        'proj=longlat'
        'ellps=WGS84'
        'towgs84=0,0,0,0,0,0,0'
        'no_defs'
    END
END

Best Answer

You are doing it right and 1 is the status code (SUCCESS). Below is a complete (working) example of queryByShape implemented in Java mapscript (for a vector layer).

Translating the sample into python should be a piece of cake, but I'm not sure how it works with raster layers.

/**
 * query a single layer by specifying a shape (in map coordinates).
 *
 * @param i              the layer to quer
 * @param s the shape whithin which results must be found
 * @return the query results
 */
protected List queryLayerByShape(int i, shapeObj s) {
    List results = new ArrayList();

    layerObj layer = map.getLayer(i);
    if (layer != null && isVisible(layer)) {
        if (log.isDebugEnabled()) {
            log.debug("Querying (by shape) layer: " + layer.getName());
        }
        if (layer.queryByShape(map, s) == mapscript.MS_SUCCESS) {
            if (layer.open() == mapscript.MS_SUCCESS) {
                for (int j = 0; j < layer.getNumResults(); j++) {
                    if (log.isDebugEnabled())
                        log.debug("Layer " + layer.getName() + ", result number:" + j);
                    resultCacheMemberObj resultMember = layer.getResult(j);
                    shapeObj shape = new shapeObj(layer.getType());
                    layer.getShape(shape, resultMember.getTileindex(), resultMember.getShapeindex());
                    if (shape != null) {
                        Map aResult = getResult(shape, layer);

                        results.add(aResult);
                    } else {
                        log.error("Shape " + j + " is null!");
                    }
                }
        // not necessary since 5.6
                //layer.close();
            } else {
                log.error("Cannot open layer: " + layer.getName());
            }
        } else {
            log.info("Query by shape on layer " + layer.getName() + " failed.");
        }
    }
    return results;
}

public shapeObj getSelectionShape(MapRequest request) {
    shapeObj s=new shapeObj(mapscript.MS_SHAPE_POLYGON);
    double x,y;

    String[] coords=request.getSelectionCoordinates().split(";");
    lineObj l=new lineObj();
    for ( int i=0; i<coords.length; i++ ) {
        String[] pair=coords[i].split(",");
        x=XToMapCoordinates(map.getWidth(), map.getExtent(), Double.parseDouble(pair[0]));
        y=YToMapCoordinates(map.getHeight(), map.getExtent(), Double.parseDouble(pair[1]));

        l.add(new pointObj(x, y, 0));
    }
    s.add(l);

    if (log.isTraceEnabled()) {
        log.trace("QueryByShape selectionCoords="+request.getSelectionCoordinates());
    }

    return s;
}

Mapserver documentation about raster querying:

http://trac.osgeo.org/mapserver/wiki/RasterQuery

http://mapserver.org/input/raster.html#raster-query

Related Question