1

I am currently working on code that will determine the connectivity of pixels within a raster image. I am currently having an issue with selecting all pixels within a minimum distance of a given pixel.

The goal of this line of code is to extract a subset from a matrix whose entries look like this: 

             jcol      icol    valcol
300   0.005076142 0.8860759 0.0000000
27077 0.862944162 0.2911392 0.0000000
30604 0.974619289 0.4746835 0.0000000
3138  0.096446701 0.7341772 0.0000000
9926  0.314720812 0.4240506 0.0000000
27328 0.868020305 0.8734177 0.0000000
17624 0.558375635 0.8417722 0.0000000
13117 0.416243655 0.4936709 0.0000000
22622 0.720812183 0.2721519 0.6509804
15720 0.497461929 0.8670886 0.0000000

This where jcol is the relative y position, icol is the relative x position and valcol is the value of the raster at the pixel (x,y). This matrix covers the entire raster.

Currently I am attempting to use the subset function to go about this task:

P_N = subset(cluster.matrix, abs(cluster.matrix[,1]-x_N[,1]) == xthresh | abs(cluster.matrix[,2]-x_N[,2]) == ythresh)

Where xhtresh and ythresh represent the distance around the sample point x_N, a point which is chosen from the original (example below):

     jcol      icol    valcol 
3099    1 0.6518987 0.4862745

The issue is that I get the following output for P_N:

[1] jcol     icol     valcol
<0 rows> (or 0-length row.names)

Ideally, P_N would return the 8 points directly adjacent to x_N

I think my syntax in the subset function is wrong, especially because x_N is a 1x4 data entry, while cluster.matrix is a ~7000x4 matrix.

JasonAizkalns
  • 20,243
  • 8
  • 57
  • 116
TimColeman
  • 21
  • 3
  • 1
    Might `raster::clump` [(see e.g. here)](http://stackoverflow.com/questions/20659186/combining-polygons-and-calculating-their-area-i-e-number-of-cells-in-r/20663885#20663885) help with your larger question, about connectivity of pixels in the raster? – Josh O'Brien Feb 23 '16 at 19:02
  • I suggest that you replace the "==" signs with "<" and try again. – RHertel Feb 23 '16 at 19:05
  • 1
    On the lines of `clump`, see also [this gist](https://gist.github.com/johnbaums/6d155cfca28e02b05ad5), for which you can specify a neighbourhood distance (ping @JoshO'Brien). – jbaums Feb 23 '16 at 20:40
  • @jbaums -- Looks intriguing, thanks. Unfortunately, I can't so far get it to work on my Windows box, even though I have found [this alternate version](https://gist.github.com/Pakillo/c6b076eceb0ef5a70e3b) of `polygonizer.R` :(. Will come back to this, though, when I have some time. – Josh O'Brien Feb 24 '16 at 04:05
  • @JoshO'Brien - yep, the GDAL/Python/R interaction can be a bit fussy... :/ – jbaums Feb 24 '16 at 04:39
  • 1
    @jbaums I found a fix, and made the edit in [my own fork of that gist](https://gist.github.com/JoshOBrien/2bf028356f38538b5ddf). (Basically avoid the errors I was getting, I needed to pretty up `pypath` in a couple of ways, by adding this line: `pypath <- tools::file_path_sans_ext(normalizePath(pypath))`. I'm not real familiar with how gists work, though. Is there some way you can incorporate that fix and Pakillo's in your original gist? – Josh O'Brien Feb 24 '16 at 17:34
  • nice work, @JoshO'Brien - thank you. I'll merge the changes. Did you need to use OSGeo4W? – jbaums Feb 24 '16 at 20:28
  • 1
    @jbaums Hi John. Yep, that's what I used. (I debugged `polygonizer()` by first opening a OSGeo4W and figuring out what command it wanted and then worked back from there. And in fact, on my machine at least, simply doing `pypath <- "gdal_polygonize"` works just fine. I'd kind of expect it to also do so on other Windows machines that have OSGeo4W64 installed, though there doesn't seem to be much harm in having it be looked for on the search path by `Sys.which()`.) Thanks again for pointing me to those gists! – Josh O'Brien Feb 24 '16 at 20:45

1 Answers1

2

I think you want raster::adjacent.

library(raster)
r <- raster(matrix(runif(100), 10))
cells <- adjacent(r, 34, 8, pairs=FALSE)
adj <- cbind(xyFromCell(r, cells), value=extract(r, cells))

Above, adjacent returns the cell numbers of the 8 cells adjacent to cell 34 of raster r. You can also pass a neighbourhood matrix that you might derive from a distance matrix/raster - see the help for argument directions at ?adjacent.

We then extract the x and y coordinates using xyFromCell, and the cell values with extract (alternatively we could just use r[cells]).

jbaums
  • 27,115
  • 5
  • 79
  • 119
  • Thanks for the suggestion -- when using `adjacent` how does one identify a particular cell using the `cells` argument? The convention in the description seems vague. – TimColeman Feb 24 '16 at 02:28
  • @TimColeman - Are you asking how to specify the cells you want to check adjacency for? You can use `cellFromXY` if you have coords and need to get the cell number. – jbaums Feb 24 '16 at 03:03