Question re: PostGIS and MapScript getExtent() function

[Steve Hall]

Hi Paul,

Exactly right - it definitely needs reviewing by those that know more about 
MapServer internals than I do.  (which is most of you I expect).

How is this function used my MapServer itself then?  The context in which I 
am using it is via an explicit $layerObj->getExtent() call from 
PHP/MapScript, so in this respect the current implementation is not much 
use, but I suspect there is more to this function than just an 
implementation for this call.   The main issue I had was that this function 
worked as I expected for Oracle, yet not for PostGIS, meaning I'd need some 
ugly code in my app to deal with the inconsistency...

Anyway, i've attached my version of the function for you to take a look at. 
My implemenation is basically the version that was commented out, with the 
compilcation errors & warning resolved, added support for the FILTER 
expression and support for either BOX or BOX3D return clauses from PostGIS 
extent().

Many thanks,

Steve



----- Original Message ----- 
From: "Paul Ramsey" <pramsey at refractions.net>
To: "Steve Hall" <steve.hall at MERCATORGEOSYSTEMS.CO.UK>
Cc: <MAPSERVER-USERS at LISTS.UMN.EDU>
Sent: Tuesday, July 12, 2005 3:49 PM
Subject: Re: [UMN_MAPSERVER-USERS] Question re: PostGIS and MapScript 
getExtent() function


> msPOSTGISLayerGetExtent() is potentially extremely inefficient,  depending 
> on how it is implemented, so I would have to see your patch  to judge if 
> it does more harm than good. It may be doing good things  for you and bad 
> things for others.  One thing it *could* do without  harm is read the 
> wms_extent meta-parameter and use that if it is  provided by the user. 
> There are also a few version specific fast  extent functions available in 
> PostGIS, but they are not guaranteed to  work in all cases, at all times, 
> correctly.
>
> On 12-Jul-05, at 2:58 AM, Steve Hall wrote:
>
>> Hi All,
>>
>> I'm using MapServer 4.6.0 with Oracle Spatial and PostGIS data  sources 
>> via PHP/MapScript.  Part of my application code looks like  this 
>> pseudo-code :
>>
>>     $map->setLayerFilter("ThematicData", "REQUEST_ID=121");
>>     $map->zoomToLayer("ThematicData");
>> Here i'm setting a specific FILTER expression for a layer, and then 
>> zooming to the extent of the (now filtered) layer.  I use the 
>> layerObj->getExtent() and then mapObj->setExtent() functions to  achieve 
>> this behind the scenes.  This works great for Oracle, the  query is fast 
>> and makes for nice, tidy code in my application  without the need for 
>> Oracle specific SQL statements to derive the  new extent.
>>
>> However when I try this on PostGIS I just get a massive (useless)  extent 
>> returned.  Digging into the code I can see that the 
>> msPOSTGISLayerGetExtent() function in mappostgis.c is effectively  doing 
>> nothing although a lot of code is commented out.
>>
>> Working on my local code I have re-enabled this code, fixed the 
>> compilation errors and changed it so that when a FILTER is applied  to 
>> the layer this is used as part of the SQL call to PG.  It seems  to work 
>> great (from PHP/MapScript) and performance is fine and is  really useful 
>> for me meaning I don't need to calculate the extents  for a map myself - 
>> MapServer can supply the information in a  generic and abstract fashion.
>>
>> Ideally I would like to get this code fed back into MapServer, but  of 
>> course only if it's acceptable to you all.  So,  should I  instead post 
>> this to mapserver-dev or can I pass on the new code to  someone from 
>> here?
>>
>> There may be good, historical, reasons why this wasn't ever  implemented 
>> though so please let me know!
>>
>> Many thanks,
>>
>> Steve Hall
>>
>> 
-------------- next part --------------
/* This code executes the PostGIS extent() function on the supplied DATA statement */
/* and filter (if supplied).  */
int msPOSTGISLayerGetExtent(layerObj *layer, rectObj *extent)
{
    if(layer->debug) {
        msDebug("msPOSTGISLayerGetExtent called\n");
    }

    PGresult   *query_result;
    char        sql[5000];

    msPOSTGISLayerInfo *layerinfo;

    char    * table_name       = 0;
    char    * geom_column_name = 0;
    char    * urid_name        = 0;
    char    * user_srid        = 0;
    if (layer == NULL)
    {
        char tmp[5000];
        sprintf(tmp, "layer is null - have you opened the layer yet?");
        msSetError(MS_QUERYERR, tmp, "msPOSTGISLayerGetExtent()");

        return(MS_FAILURE);
    }

    /*  get query components */
    layerinfo = (msPOSTGISLayerInfo *) layer->layerinfo;
    msPOSTGISLayerParseData(layer, &geom_column_name, &table_name, &urid_name, &user_srid, layer->debug);
    /* append filter string, if supplied, so that we're only getting the extent of the required data */
    if (layer->filter.string != NULL) {
        sprintf(sql, "select extent(%s) from %s where %s", geom_column_name, table_name, (layer->filter.string));
    } else {
        sprintf(sql, "select extent(%s) from %s", geom_column_name, table_name);
    }
    if(layer->debug) {
        msDebug("msPOSTGISLayerGetExtent SQL statement is %s\n", sql);
    }

    /* check for a connection */
    if (layerinfo->conn == NULL)
    {
        char tmp[5000];
        sprintf(tmp, "layer doesnt have a postgis connection - have you opened the layer yet?");
        msSetError(MS_QUERYERR, tmp, "msPOSTGISLayerGetExtent()");

        return(MS_FAILURE);
    }

    /* Execute the query */
    query_result = PQexec(layerinfo->conn, sql);
    if (!(query_result) || PQresultStatus(query_result) !=  PGRES_TUPLES_OK)
    {
        char tmp[5000];
        sprintf(tmp, "Error executing POSTGIS  SQL   statement (in msPOSTGISLayerGetExtent): %s", layerinfo->sql);
        msSetError(MS_QUERYERR, tmp, "msPOSTGISLayerGetExtent()");

        PQclear(query_result);
        return(MS_FAILURE);
    }

    /* Check only one row was returned */
    if (PQntuples(query_result) != 1)
    {
        char tmp[5000];
        sprintf(tmp, "Error executing POSTGIS  SQL   statement (in msPOSTGISLayerGetExtent) [doesnt have exactly 1 result]: %s", layerinfo->sql);
        msSetError(MS_QUERYERR, tmp, "msPOSTGISLayerGetExtent()");

        PQclear(query_result);
        return(MS_FAILURE);
    }

    /* Parse out extent from returned "BOX" (or BOX3D) clause */
    char * qry_extent = PQgetvalue(query_result,0,0);
    if (strncmp(qry_extent, "BOX3D(", 6) == 0) {
		    // Z axis values are retrieved but discarded
		    double minZ, maxZ;
	      sscanf(qry_extent, "BOX3D(%lf %lf %lf,%lf %lf %lf)", &extent->minx, &extent->miny, &minZ, &extent->maxx, &extent->maxy, &maxZ );
	  } else if (strncmp(qry_extent, "BOX2D(", 6) == 0) {
  		  sscanf(qry_extent, "BOX2D(%lf %lf,%lf %lf)", &extent->minx, &extent->miny, &extent->maxx, &extent->maxy );
	  } else if (strncmp(qry_extent, "BOX(", 4) == 0) {
		    sscanf(qry_extent, "BOX(%lf %lf,%lf %lf)", &extent->minx, &extent->miny, &extent->maxx, &extent->maxy );
	  } else {
        char tmp[5000];
        sprintf(tmp, "Error parsing returned POSTGIS extent %s", qry_extent);
        msSetError(MS_QUERYERR, tmp, "msPOSTGISLayerGetExtent()");

        PQclear(query_result);
        return(MS_FAILURE);
	  }

    PQclear(query_result);
    return MS_SUCCESS;
}