[postgis-devel] Spatial relationship functions

[Bborie Park]

On Mon, Jul 23, 2012 at 6:46 AM, Pierre Racine
<Pierre.Racine at sbf.ulaval.ca> wrote:
> Bborie,
>
>> You've got a valid point.  I'll revert based upon the ongoing thread
>> about the future of raster/geometry and geometry/raster spatial
>> relationship function variants?
>
> Did you do some tests? This implementation is faster when there is not so many nodata values but very slow when there are a lot. This is why I planned a IsNodataValue flag meaning "this raster is ONLY nodata values". If this flag is set ST_Intersects can return false right away. I think in many cases most tiles containing at least one no data value are just nodata value. They are generally tiles filling the area outside the area of interest of the data set. So I think it is a good idea to do such an optimization. However this flag is hardly updated when appropriate... Still lots of work to do to make it work properly.
>

No, I haven't tried any testing.  The two scenarios you mention are on
opposite extremes, particularly if the tiles are big.  My guess would
be that for a specific set of tiles, GDAL's polygonize will minimize
the gap between the extremes.

>> I'm currently leaning towards Paul's point of view on eliminating the
>> raster/geometry and geometry/raster spatial relationship function
>> variants...
>>
>> 1. We don't need raster-space variants.  Converting a geometry to a
>> raster loses accuracy vs converting the raster to a geometry with
>> ST_Polygon.
>>
>> 2. Besides of the performance gains from ST_Intersects(geometry,
>> raster) using your implementation, none of the other geometry/raster
>> spatial relationship functions can operate in a similar manner.
>
> If this is really true then I have no problem eliminating/not implementing them.
>

I'm certain that this is true.  For item #1, this is absolutely true
for any geometry that isn't a rectangle, say a triangle or a circle.
Rasterizing a triangle or a circle results in a raster whose area does
not match the input geometry.  Vectorizing a raster for comparison
though ensures that the geometry of the raster's band is spatially
correct.

#2 is true because there really is no clean mathematical way for any
of the spatial relationships besides ST_Intersects.  Therefore, we
have to operate in the geometry world when testing for relationships.

> Did you try, to optimize ST_Polygon, to set all the withvalue pixels to a unique value before passing the raster to GDAL for vectorization? I'm pretty sure the GDAL vectorizer is much faster when the raster is more homogenous.
>

The problem with setting all the non NODATA value pixels to a uniform
value is that you'd have to iterate over all the pixels.  GDAL's
polygonize (both int and float variants) also have to iterate over the
pixels... so no real gain here.

>> Calling ST_SomeTest(ST_Polygon(raster), geometry) is what any
>> geometry/raster spatial relationship function would essentially be.
>> Less functions = less maintenance.  And with Regina indicating the
>> same issues with geometry/geography and how they addressed that (make
>> users explicitly convert the geography to a geometry), I'm fine with
>> doing the same with geometry/raster.
>
> But we keep all the raster/raster ones right? And we keep ST_Intersects(geometry/raster)?
>

Correct.  ST_Intersects(raster, raster), ST_Intersects(geometry,
raster) and ST_Intersects(raster, geometry) would both remain.

> How/where shall we document that ST_SomeTest(ST_Polygon(raster), geometry) is the way to do?
>

In a "Notes" section for each function.

> I have to tell you that since my boss does not have any more money to pay my work on PostGIS I had to find something else. I'm glad I could find another project here for the coming year working on agent based modeling (which I find way more challenging/interesting than working on boring rasters ;-)). I will continue to contribute to PostGIS raster mainly by advising students (like the GSoC one I have right now and another one who is going to work on interpolation this fall) and writing in my blog (coming soon). As a last task I will update the specifications pages.
>

Sorry to hear that funding ran out.  I totally understand the hunt for
more money :-).  I do expect to be working on PostGIS Raster for the
foreseeable future as the next major thing I needs of PostGIS Raster
is modelling with MapAlgebra on N-rasters.

One thing I did notice when I was digging into what other
raster-processing software (ArcGIS, ENVI, ERDAS Imagine) provides for
spatial relationships is that they don't provide any.  Well, any that
returns a boolean.

-bborie


-- 
Bborie Park
Programmer
Center for Vectorborne Diseases
UC Davis
530-752-8380
bkpark at ucdavis.edu