Last active
September 9, 2019 03:01
-
-
Save johnbaums/6d155cfca28e02b05ad5 to your computer and use it in GitHub Desktop.
Clump raster cells into patches, with optional neighbourhood distance.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
patchify <- function(x, distance, p4s, givedist=TRUE) { | |
# x: a binary Raster layer (0 or NA for background, and 1 for areas to be clumped) | |
# distance: the neighbourhood distance. Patches that occur within this distance of | |
# one another will be clumped. This should be in the units of the CRS given | |
# in p4s. If this is zero, the function will identify patches defined by | |
# 8-connectivity (i.e. queen's case). | |
# p4s: an equal-area projection appropriate for the region. This will be used when | |
# calculating inter-patch distances. The returned objects will be in the | |
# original coordinate system. | |
# givedist: should the distance matrix be returned? (logical). Distances are in the | |
# units of p4s, and are shortest cartesian distances between patches. | |
require(raster) | |
require(rgeos) | |
require(SDMTools) | |
if(!is(x, 'Raster')) x <- raster(x) | |
if(!is(p4s, 'CRS')) p4s <- CRS(p4s) | |
if(is.na(proj4string(x))) stop(substitute(x), ' lacks a CRS.') | |
x[x == 0] <- NA | |
cc <- ConnCompLabel(x) | |
p <- polygonizer(cc) # https://gist.github.com/johnbaums/26e8091f082f2b3dd279 | |
suppressWarnings(proj4string(p) <- proj4string(x)) | |
pproj <- spTransform(p, p4s) | |
bproj <- gBuffer(pproj, width=distance) | |
pproj$patch <- over(pproj, disaggregate(bproj)) | |
b <- spTransform(pproj, CRS(proj4string(x))) | |
pout <- aggregate(b, by='patch') | |
pout$patch <- as.factor(pout$patch) | |
rout <- rasterize(pout, x) | |
out <- list(patchrast=rout, patchpoly=pout) | |
if(isTRUE(givedist)) { | |
poutproj <- spTransform(pout, p4s) | |
d <- gDistance(poutproj, poutproj, byid=TRUE) | |
out <- c(out, list(distance=d)) | |
} | |
return(out) | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Example:
First, simulate some spatially-correlated data
Source the required gists
Here we load the
polygonizer
(see here) andpatchify
gists.Now identify patches
Here we set the neighbourhood distance to 1000 m, and specify Australian Albers (i.e., EPSG:3577) as the projection that will be used when calculating inter-patch distances. (NB: The returned raster and polygon objects will revert to the original coordinate system.)
The function returns the patches as a
Raster
object and as aSpatialPolygonsDataFrame
. Below, we plot each.What happens if we increase the neighbourhood distance to 10 kilometres?