Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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
share|improve this question
add comment

1 Answer

up vote 3 down vote accepted

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

share|improve this answer
    
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. –  Calvin Cheng Dec 27 '11 at 8:39
    
do you call layer.open AND have a template in your layer definition? Without a template querying will not work. –  unicoletti Dec 27 '11 at 9:22
    
I have a template defined as follows:- In [14]: layer.template Out[14]: 'fooOnlyForWMSGetFeatureInfo' –  Calvin Cheng Dec 28 '11 at 4:01
    
Added the LAYER definition in my mapfile in my original question above. –  Calvin Cheng Dec 28 '11 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 Dec 28 '11 at 6:19
show 4 more comments

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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