raster

r = raster(matrix(sample(1:9, 100, replace = TRUE), 10, 10))

write

writeRaster(r, "raster.tif")

read

r = raster("raster.tif")

Access

Direct index r[i]

Using line and columns r[lin, col]

Add vectors

cellFromXY(r, pt)

## [1] 23

colFromX(r, pt$x)

## [1] 3

rowFromY(r, pt$y)

## [1] 3

fourCellsFromXY(r, as.matrix(pt))

##      [,1] [,2] [,3] [,4]
## [1,]   23   33   32   22

cellFromLine(r, ln)

## [[1]]
##  [1] 35 36 45 55 61 62 63 64 65 71 81 91

cellFromPolygon(r, poly)

## [[1]]
##  [1] 24 25 26 34 35 36 37 45 46 47 55 56 57 65 66 67 75 76 77

aggregate

r_agg = aggregate(r, 2)

disaggregate

r_dis = disaggregate(r, 2)

cover

r_missing = raster(matrix(NA, 10, 10))
r_missing[41:60] = 6

r_covered = cover(r_missing, r)

mask

r_masked = mask(r, r_missing)

calc

fun = function(x) { x * 10 }
r_mul = calc(r, fun)

overlay

r1 = init(r, function(x) sample(5:6, x, replace = TRUE))
r2 = init(r, function(x) sample(1:2, x, replace = TRUE))
r3 = overlay(r1, r2, fun=function(x,y){return(x+y)})

distance

We need a raster containing origins to compute distances (so we need units)

r_origin = raster(matrix(NA, 10, 10))
extent(r_origin) = extent(c(0, 10, 0, 10))
projection(r_origin) = "+proj=tmerc +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +units=m"
r_origin[c(15, 16, 45)] = 1

r_dist = distance(r_origin) 

reclassify

m = cbind(
  from    = c(0, 2, 4), 
  to      = c(2, 4, 6), 
  becomes = c(1, 2, 3)
  )

r_reclass = reclassify(r_dist, m)

focal

r_focal = focal(r, fun = mean, w = matrix(1/9,nrow=3,ncol=3))

You can remove edge effect by ignoring NA values

r_focal_nedge = focal(r, fun = mean, w = matrix(1/9,nrow=3,ncol=3), na.rm = TRUE, pad = TRUE)

zonal

zonal(r, z, mean)

clump

r_clumps = clump(r_bin)

boundaries

require("igraph")
r_clus = init(r, function(x) ifelse(runif(x) > 0.2, 1, NA))
r_boundaries = raster::boundaries(r_clus)

extract

extract(r, pt)

##   
## 2

extract(r, ln)

## [[1]]
##  [1] 3 7 4 3 6 8 1 2 1 4 6 5

extract(r, poly)

## [[1]]
##  [1] 2 7 7 2 3 7 9 4 5 8 3 8 7 1 7 1 9 5 3

rasterize

rasterize(pt, r)
rasterize(ln, r)
rasterize(poly, r)

Note that polygon rasterization is by default looking at cell centroid overlap,

density

density(r)

histogram

hist(r)

Spatial autocorrelation

Moran(r)
Geary(r)

extent

extent(r)

crop

f = extent(0.23, 0.86, 0.22, 0.73)

crop(r, f)

intersect

g = extent(0.43, 1.1, -0.1, 0.53)

intersect(r, g)

cellsFromExtent(r, g)

##  [1]  55  56  57  58  59  60  65  66  67  68  69  70  75  76  77  78  79
## [18]  80  85  86  87  88  89  90  95  96  97  98  99 100

union

union(extent(r), g)

flip

flip(r_split, "y")

extend

extend(r, 3)

mosaic

mosaic(r, r1, r2, fun = sum)

projection

projection(r)

rotate

rotate(r)

rasterToPoints

pts = rasterToPoints(r)

rasterToPolygons

pols = rasterToPolygons(r_clumps, dissolve = TRUE)

terrain

r_terrain = terrain(r_volcano, opt = c("slope", "aspect", "tpi", "tri", "roughness", "flowdir"))

Angular data can sometimes be better expressed using a circular palette. In the following figure, blue is North orientation, while South is red. Both colors reach black at West and white at East. The preceeding figure had some sharp edges on North faces, when angle slightly changed from 360 to 0.

hillShade

r_hill = hillShade(r_terrain[["slope"]], r_terrain[["aspect"]])

interpolate

tps = Tps(survey[, c("x", "y")], survey$layer)

## Warning: 
## Grid searches over lambda (nugget and sill variances) with  minima at the endpoints: 
##   (GCV) Generalized Cross-Validation 
##    minimum at  right endpoint  lambda  =  1.406207e-05 (eff. df= 47.5 )

r_interpolate = interpolate(r_volcano, tps)