[GIS] How to extract values from a raster using R

cartographyrrasterstretch

I'm a beginner at using R software, and I'm attempting to perform a percent linear stretch on the following raster data:

> KBnew

#class: RasterLayer
#dimensions  : 27018, 26389, 712978002  (nrow, ncol, ncell)
#resolution  : 15, 15  (x, y) 
#extent  : 541792.5, 937627.5, 5454772, 5860042 (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0  
#data source : /private/var/folders/sy/qh_q0tbs29l_j9wfndb6cw5r0000gp/T/RtmpDc7N7c/raster/r_tmp_2016-02-15_212425_11368_77493.grd
#names       : layer  
#values      : 7.712379, 205.2743  (min, max)

… although (if I'm doing this correctly) from what I understand, I need to know the mean and standard deviation of the raster values (which would be contained within the range of 7.712379,205.2743 or min/max values). How do I extract these values, within the minimum/maximum range?

Best Answer

The mean and standard deviation can be calculated using cellStats() from the raster package.

#install.packages(raster)
library(raster)

# example dataset
r <- raster(nrow=18, ncol=36)
r[] <- runif(ncell(r)) * 10

# get mean
cellStats(r,'mean')

Load libraries and example data:

# Load libraries
library('raster')
library('rgdal')

# Load a SpatialPolygonsDataFrame example
# Load Brazil administrative level 2 shapefile
BRA_adm2 <- raster::getData(country = "BRA", level = 2)

# Convert NAMES level 2 to factor 
BRA_adm2$NAME_2 <- as.factor(BRA_adm2$NAME_2)

# Plot BRA_adm2
plot(BRA_adm2)
box()

# Define RasterLayer object
r.raster <- raster()

# Define raster extent
extent(r.raster) <- extent(BRA_adm2)

# Define pixel size
res(r.raster) <- 0.1

Figure 1: Brazil SpatialPolygonsDataFrame plot

Simple thread example

# Simple thread -----------------------------------------------------------

# Rasterize
system.time(BRA_adm2.r <- rasterize(BRA_adm2, r.raster, 'NAME_2'))

Time in my laptop:

# Output:
# user  system elapsed 
# 23.883    0.010   23.891

Multithread thread example

# Multithread -------------------------------------------------------------

# Load 'parallel' package for support Parallel computation in R
library('parallel')

# Calculate the number of cores
no_cores <- detectCores() - 1

# Number of polygons features in SPDF
features <- 1:nrow(BRA_adm2[,])

# Split features in n parts
n <- 50
parts <- split(features, cut(features, n))

# Initiate cluster (after loading all the necessary object to R environment: BRA_adm2, parts, r.raster, n)
cl <- makeCluster(no_cores, type = "FORK")
print(cl)

# Parallelize rasterize function
system.time(rParts <- parLapply(cl = cl, X = 1:n, fun = function(x) rasterize(BRA_adm2[parts[[x]],], r.raster, 'NAME_2')))

# Finish
stopCluster(cl)

# Merge all raster parts
rMerge <- do.call(merge, rParts)

# Plot raster
plot(rMerge)

Figure 2: Brazil Raster plot

Time in my laptop:

# Output:
# user  system elapsed 
# 0.203   0.033   8.688

More info about parallelization in R:

Best Answer

Related Solutions

[GIS] How to read *.adf files into R

Process Vector to Raster Faster with R – Speed Optimization Techniques

Load libraries and example data:

Simple thread example

Multithread thread example

Related Question