[postgis-users] Raster Points in Polygon (was Re: PostGIS Raster for Met Data)

[Michael Akinde]

Hi,

I thought to provide a short follow-up to this discussion to explain how we ended up solving this issue.

We experimented a good deal with various solutions. Pierre's solution was the template that we worked the most with, adding the additional tweak of taking the bounding box of the polygon first, and then using the bounding box to limit the number of points considered in the ST_Intersects. For most polygons, this increased performance of the query by an order of magnitude or better.

We did run into a couple of issues:

- As mentioned earlier in the thread, we work with multiple SRIDs in the same table. Since new SRIDs can be added and removed more or less dynamically in the database, splitting this up into multiple tables is basically not an option. The idea of modifying PostGIS to allow multiple SRIDs in a table did not seem very appealing to us.

- If I understand correctly, PG Raster seems to be locked into a left-upper X-Y coordinate system (image), whereas we prefer the scientific standard left-lower X-Y origin system. Not a deal breaker, but it adds a layer of additional work for our conversions.

Also, performance is always a concern. The solution we have settled on for the time being is based on stored procedures built in C using GEOS. Our functions do essentially the same stuff as proposed by Pierre (though we transform the raster points to the polygon's SRID rather than the other way around), with the added tweak of the bounding box. The bounding box can - in the worst case - result in some points being omitted from the intersects which should have been included, but it's a trade off that we are willing to accept in return for performance.

With the C version, we were able to take the performance of the point-in-polygon algorithm itself from the ~40s of the brute-force version down to ~50ms (many times faster than what we could manage with plpgsql. Our benchmark returns ~7K data values/second even before we spent any significant time optimizing it, which makes us very happy, although it does mean that we probably won't integrate with Postgis raster at this time.

Sadly, our code is probably too application specific to be of use to anyone else, but it's all available on our public github if anyone is interested (https://github.com/wdb/wdb/tree/master/src/callInterface/core/wciRead/extractGridData).

Regards,

Michael A.
Database Architect, Met.no

----- Original Message -----
> Hi,
> 
> ----- Original Message -----
> > I had to modify the coordinates of your raster and your polygon
> > because the 900917 SRID do not exist by default in spatial_ref_sys,
> > but this example returns the x and y coordinates (a well as the
> > point
> > geometry and the value) of every pixels intersecting a polygon.
> 
> Apologies for the hassle; it's of course a local SRID of the kind used
> by our application (in this case +proj=ob_tran +o_proj=longlat
> +lon_0=-24 +o_lat_p=23.5 +a=6367470.0 +no_defs').
> 
> > I had to create the ST_PixelAsPoints(rast raster, band integer)
> > function which is a derivate of the ST_PixelAsPolygons(rast raster,
> > band integer) function available in this page:
> > http://trac.osgeo.org/postgis/wiki/WKTRasterUsefulFunctions
> 
> Many, many thanks - this is pretty much exactly what I was looking
> for. :-)
> 
> I had hoped that there was a better "native" way of getting the points
> inside the polygon than to output the points and running an
> ST_Intersects over them. Hmm - well, at least it is easy enough to
> optimize ST_PixelAsPoints by constraining it to only return the points
> within the axis-aligned minimum bounding box of the polygon. Also, I
> expect simply constructing geom using the corner point and dX/dY would
> be a bit faster than using ST_Centroid(ST_PixelAsPolygon), so perhaps
> we can get a solution along these lines running at a reasonable clip.
> 
> Thanks again for pointing us at the referenced functions - this gives
> us something to work with.
> 
> Regards,
> 
> Michael A.
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users