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

1 Answer 1

3

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

9
  • I have a python equivalent of what you wrote above in the for loop. However, I seem to be having trouble with layer.getNumResults() giving me a 0. Commented Dec 27, 2011 at 8:39
  • do you call layer.open AND have a template in your layer definition? Without a template querying will not work.
    – unicoletti
    Commented Dec 27, 2011 at 9:22
  • I have a template defined as follows:- In [14]: layer.template Out[14]: 'fooOnlyForWMSGetFeatureInfo' Commented Dec 28, 2011 at 4:01
  • Added the LAYER definition in my mapfile in my original question above. Commented Dec 28, 2011 at 4:11
  • It just occurred to me that the layer STATUS is OFF. Do you change it in Python mapscript to DEFAULT or ON before querying? Query engine skips off layers. Try changin in the mapfile STATUS=DEFAULT
    – unicoletti
    Commented Dec 28, 2011 at 6:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.