Patchify and Polygonize for Linux (when python is installed using anaconda)

patchify_polygonize_linux.r

##Not my functions but slighly mod for my use

patchify2<-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 <- polygonizer2(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)
}



polygonizer2<-function(x, outshape=NULL, pypath=NULL, readpoly=TRUE, 
         fillholes=FALSE, aggregate=FALSE, 
         quietish=TRUE) {
  # x: an R Raster layer, or the file path to a raster file recognised by GDAL 
  # outshape: the path to the output shapefile (if NULL, a temporary file will 
  #           be created) 
  # pypath: the path to gdal_polygonize.py or OSGeo4W.bat (if NULL, the function 
  #         will attempt to determine the location)
  # readpoly: should the polygon shapefile be read back into R, and returned by
  #           this function? (logical) 
  # fillholes: should holes be deleted (i.e., their area added to the containing
  #            polygon)
  # aggregate: should polygons be aggregated by their associated raster value?
  # quietish: should (some) messages be suppressed? (logical)
  if (isTRUE(readpoly) || isTRUE(fillholes)) require(rgdal)
  if (isTRUE(aggregate)) require(rgeos)
  if (is.null(pypath)) {
    cmd <- Sys.which('OSGeo4W.bat')
    pypath <- 'gdal_polygonize'
    if(cmd=='') {
      cmd <- '/home/wienslab/anaconda3/bin/python' ##This needs to be changed when python is installed using anaconda
      #pypath <- Sys.which('gdal_polygonize.py')
      pypath <- '/home/wienslab/anaconda3/bin/gdal_polygonize.py' ##This needs to be changed when python is installed using anaconda
      
      if (!file.exists(pypath)) 
        stop("Could not find gdal_polygonize.py or OSGeo4W on your system.") 
    }
  }
  if (!is.null(outshape)) {
    outshape <- sub('\\.shp$', '', outshape)
    f.exists <- file.exists(paste(outshape, c('shp', 'shx', 'dbf'), sep='.'))
    if (any(f.exists)) 
      stop(sprintf('File already exists: %s', 
                   toString(paste(outshape, c('shp', 'shx', 'dbf'), 
                                  sep='.')[f.exists])), call.=FALSE)
  } else outshape <- tempfile()
  if (is(x, 'Raster')) {
    require(raster)
    writeRaster(x, {f <- tempfile(fileext='.tif')})
    rastpath <- normalizePath(f)
  } else if (is.character(x)) {
    rastpath <- normalizePath(x)
  } else stop('x must be a file path (character string), or a Raster object.')
  
  system2(cmd, args=(
    sprintf('"%s" "%s" %s -f "ESRI Shapefile" "%s.shp"', 
            pypath, rastpath, ifelse(quietish, '-q ', ''), outshape)))
  
  if(isTRUE(aggregate)||isTRUE(readpoly)||isTRUE(fillholes)) {
    shp <- readOGR(dirname(outshape), layer=basename(outshape), 
                   verbose=!quietish)    
  } else return(NULL)
  
  if (isTRUE(fillholes)) {
    poly_noholes <- lapply(shp@polygons, function(x) {
      Filter(function(p) p@ringDir==1, x@Polygons)[[1]]
    })
    pp <- SpatialPolygons(mapply(function(x, id) {
      list(Polygons(list(x), ID=id))
    }, poly_noholes, row.names(shp)), proj4string=CRS(proj4string(shp)))
    shp <- SpatialPolygonsDataFrame(pp, shp@data)
    if(isTRUE(aggregate)) shp <- aggregate(shp, names(shp))
    writeOGR(shp, dirname(outshape), basename(outshape), 
             'ESRI Shapefile', overwrite=TRUE)
  }
  if(isTRUE(aggregate) & !isTRUE(fillholes)) {
    shp <- aggregate(shp, names(shp))
    writeOGR(shp, dirname(outshape), basename(outshape), 
             'ESRI Shapefile', overwrite=TRUE)
  }
  ifelse(isTRUE(readpoly), return(shp), return(NULL))
}