[raster] Two raster ST_Intersects

[Bborie Park]

Need to implement a ST_Intersects(raster, raster) upon which other ST_Intersects should call.

[bnordgren]

I don't suppose I could twist your arm into helping with the generation2 raster iterator for this, could I? I have most of the "toolbox" written, meaning the interfaces as well as "providers" for geometries and rasters. I also have "wrapper" providers, to do things like represent the object in a different projection. And don't forget the spatial and value "operator primitives". Really, the next step is to write the thing which evaluates the inputs at fixed intervals (e.g., grid cells)…Almost there…

The code is on github: ​https://github.com/bnordgren/postgis

My branches are messy because I was learning how to use git, but the relevant code is on 'ri-gen2-arch'. If you want to build it, you have to:

export PYTHON=python
./configure ...

otherwise the configure will fail.

[Bborie Park]

A proposed set of two raster ST_Intersects functions. All existing raster/geometry functions will stay the same, though they may be adjusted to decide if converting the geometry to a raster is better solution that converting the raster to a set of geometries.

1. ST_Intersects(raster rastA, integer bandA, boolean exclude_nodata_valueA, raster rastB, integer bandB, boolean exclude_nodata_valueB) → boolean

2. ST_Intersects(raster rastA, integer bandA, raster rastB, integer bandB, ) → boolean

3. ST_Intersects(raster rastA, boolean exclude_nodata_valueA, raster rastB, boolean exclude_nodata_valueB) → boolean

4. ST_Intersects(raster rastA, raster rastB, integer bandA DEFAULT 1, boolean exclude_nodata_valueA DEFAULT TRUE, integer bandB DEFAULT 1, boolean exclude_nodata_valueB DEFAULT TRUE) → boolean

5. ST_Intersects(raster rastA, raster rastB, integer bandA, integer bandB) → boolean

6. ST_Intersects(raster rastA, raster rastB, boolean exclude_nodata_valueA, boolean exclude_nodata_valueB DEFAULT TRUE) → boolean

[bnordgren]

I'm somewhat at a loss to understand why the band matters (other than to compare with the nodata value, of course). So if the user elects to "exclude_nodata_valueA", then the convex hull of the rastA is the shape which needs to be tested for intersection with rastB, right? Essentially what I'm getting at is that the only reason to supply a band index at all is to tell the function which band contains the nodata value. If you want to ignore the nodata value, I think you can omit the band index entirely.

I think this simplifies the above list to #5, #2 and

ST_Intersects(raster rastA, raster rastB) → boolean

[pracine]

I agree with Bryce on this. Could we get some more detailed specs on this operator in the wiki before going further? There is a lot of detail underlying this operator. They should be described and maybe planned in a broader context.

[Bborie Park]

I've added detailed specs of the algorithm and variations of the function in the wiki.

[bnordgren]

Oh my. That seems complicated. May I offer an alternative algorithm which should mean less operations per pixel? It may be too simple or not cover a use case, so please blast away at it. It's pseudocode so holler out if anything is unclear.

if (A.SRID != B.SRID) return false; 
grid_overlaps := convex hull of A overlaps convex hull of B.
if (no bands have been provided) { 
  return grid_overlaps; 
} else {
  if (not grid_overlaps) return false ;

  smaller := A or B, depending on which raster has fewer pixels;
  big     := the other raster
  for smaller_j := 1,small_HEIGHT {
     for smaller_i := 1, small_WIDTH { 
         (big_i, big_j) = big_geopoint2cell(small_cell2geopoint(smaller_i, smaller_j))
         if ( 0 < big_i <= big_WIDTH) AND (0 < big_j <= big_HEIGHT) {
             round big_i and big_j to nearest integer indices ;
             if (big_getValue(big_i, big_j) != NODATA) {
                return TRUE ; 
             }
         }
     }
  }          
}
return FALSE ;

This could address the special case where one raster fits completely within one cell of the other if we make sure that the selection of "smaller" does not allow the pixelsize of "smaller" to be larger than the whole "bigger" raster.

So maybe the criteria for "smaller" would be "which raster has the minimum pixel size", and if there's a tie: "which raster has fewer pixels"?

Note, I implemented the rt_raster_geopoint_to_cell(), offloaded both forward and backwards transformations to GDAL, and made sure the code uses the accessor methods instead of the fields here:

​https://github.com/bnordgren/postgis/commit/ab632411b303aca898c285c53e428060ea9e35de

[Bborie Park]

I don't believe your algorithm provides for the situation where two pixels (one from each raster) touch, such as when two rasters are next to one another. In addition, if the two rasters only intersect at one pixel, looping over every pixel to find that one pixel may take far longer than need be.

By using the grid lines of each raster in an every third line pattern and determining their intersection points, I can test for overlap, within and touches quickly and in one pass.

[bnordgren]

This can get pretty dicey pretty quick. My algorithm applies a nearest neighbor sampling of "big" using indices derived from "small". It could be that the pixel from "small" actually straddles the dividing line between two to four pixels from "big", and further, it could be "more in" the NODATA pixel and get missed. I suspect that checking adjacent pixel(s) from big which are "within" one small pixel distance could help. However, one cannot assume that the two rasters are aligned and this may throw a wrench in the works.

IIRC, using a 3x3 kernel does not typically mean you get to advance by three pixels with every iteration. Normally, you center the 3x3 kernel on every pixel, perform the operation, and move to the very next pixel. I'm not sure I understand what you're proposing to do or how it would catch intersections in the skipped pixels. We may need to start drawing pictures.

[Bborie Park]

Actually, I've already written what I've envisioned and am currently in the process of writing regression tests for it and seeing how it performs will real data. The algorithm I've written doesn't care if the rasters are aligned, just that they are of the same SRID.

I'm not using a 3x3 kernel, rather each pass uses every third grid lines. So the first pass uses grid lines 0, 3, 6, 9, …, the second 1, 4, 7, 10,… and the third 2, 5, 8, 11…

[bnordgren]

I think I'd still like to see a picture or an algorithm description. I can retire myself to poking holes in your solution rather than placing one of my own in harm's way.

[pracine]

Note, I implemented the rt_raster_geopoint_to_cell(), offloaded both forward and backwards transformations to GDAL, and made sure the code uses the accessor methods instead of the fields here:

​https://github.com/bnordgren/postgis/commit/ab632411b303aca898c285c53e428060ea9e35de

If I understand well these could be used to implement c version of all the world2raster and raster2world plpgsql functions?

[bnordgren]

Replying to pracine:

If I understand well these could be used to implement c version of all the world2raster and raster2world plpgsql functions?

Yep. It's actually performed by GDAL using the GDALApplyGeoTransform() function. GDAL is also used to calculate the inverse transform whenever the accessors are used to change the parameters of the affine transform (scale, skew, offset). In essence, the forward and reverse transforms are always pre-computed and ready for use.

[pracine]

From my point of view ST_Intersects should be a wrapper around the iterator part of the planned two raster version of mapalgebra.

For this special case the iterator should be able to traverse all the raster pixels in a pseudo sequencial way (every 0 + 3 pixels on the first pass, every 1 + 3 pixels on the second pass and every 2 + 3 pixels on the third pass (this is different than sampling where skipped pixels are forgotten forever)) since we are searching for pixel having a certain characteristic and most nodata value pixels are generally neighbors. This is the technique used by ST_Intersects(raster, geometry) and it prooved much faster than a normal scan.

For this special case, the iterator should also return as soon as both values are not nodata.

I know this seems to be a nice case for Bryce's iterator but I still need a demonstration that this iterator is more simple/flexible/fast/useful than the iterator used in the two raster ST_MapAlgebra plpgsql prototype (not unnecessarily introducing the concept of mask by the way .

[bnordgren]

Replying to pracine:

I know this seems to be a nice case for Bryce's iterator but I still need a demonstration that this iterator is more simple/flexible/fast/useful than the iterator used in the two raster ST_MapAlgebra plpgsql prototype (not unnecessarily introducing the concept of mask by the way .

The iterator part of gen2 has been written for a while, although I still need to expose it up to the SQL layer. For your application, it's main feature is that it's accessible (e.g., it's not in the middle of some other function.)

However, I was thinking about it and I don't think you want either my iterator or mapalgebra's. You want to stop iterating at the first true value, and neither iterator does that. (And I really don't think we should add that feature to either iterator.) Likewise, you don't care about storing the result for every pixel: you just need a single true/false at the end. I think these features are common to all of the DE-9M tests. This may be worth writing a second iterator which stops early and doesn't store per-pixel values.

[Bborie Park]

Bborie's implementation of two-raster ST_Intersects

[Bborie Park]

I've attached a patch that adds my implementation of ST_Intersects.

[bnordgren]

I've gotta sit down and soak this in a bit more. But I did notice a commented-out attempt to use GDALInvGeoTransform. Was there an issue with it?

[Bborie Park]

No issue with GDALInvGeoTransform. I just used it to validate the equations.

[Bborie Park]

Added in r7884