Opened 9 years ago

Closed 9 years ago

#1176 closed enhancement (fixed)

[raster] Two raster ST_Intersects

Reported by: Bborie Park Owned by: Bborie Park
Priority: medium Milestone: PostGIS 2.0.0
Component: raster Version: master
Keywords: history Cc:

Description

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

Attachments (1)

raster.diff (53.8 KB) - added by Bborie Park 9 years ago.
Bborie's implementation of two-raster ST_Intersects

Download all attachments as: .zip

Change History (20)

comment:1 Changed 9 years ago by Bborie Park

Owner: changed from pracine to Bborie Park

comment:2 Changed 9 years ago by 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.

comment:3 Changed 9 years ago by Bborie Park

Status: newassigned

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
  1. ST_Intersects(raster rastA, integer bandA, raster rastB, integer bandB, ) -> boolean
  1. ST_Intersects(raster rastA, boolean exclude_nodata_valueA, raster rastB, boolean exclude_nodata_valueB) -> boolean
  1. 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
  1. ST_Intersects(raster rastA, raster rastB, integer bandA, integer bandB) -> boolean
  1. ST_Intersects(raster rastA, raster rastB, boolean exclude_nodata_valueA, boolean exclude_nodata_valueB DEFAULT TRUE) -> boolean

comment:4 Changed 9 years ago by 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

comment:5 Changed 9 years ago by 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.

comment:6 Changed 9 years ago by Bborie Park

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

comment:7 Changed 9 years ago by 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

comment:8 Changed 9 years ago by 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.

comment:9 Changed 9 years ago by 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.

comment:10 Changed 9 years ago by 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...

comment:11 Changed 9 years ago by 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. :)

comment:12 in reply to:  7 ; Changed 9 years ago by 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?

comment:13 in reply to:  12 Changed 9 years ago by 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.

comment:14 Changed 9 years ago by 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 :-).

comment:15 in reply to:  14 Changed 9 years ago by 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.

Changed 9 years ago by Bborie Park

Attachment: raster.diff added

Bborie's implementation of two-raster ST_Intersects

comment:16 Changed 9 years ago by Bborie Park

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

comment:17 Changed 9 years ago by 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?

comment:18 Changed 9 years ago by Bborie Park

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

comment:19 Changed 9 years ago by Bborie Park

Keywords: history added
Resolution: fixed
Status: assignedclosed

Added in r7884

Note: See TracTickets for help on using tickets.