a spatial researcher's open sauce: Raster Analysis Using R

The raster package in R contains a number of interesting tools for raster analysis. It can be used within QGIS to make for easier importing of data etc. but I have been using R as standalone because of the manageR problem that necessitates a reboot. It is however not too difficult to do the work in R directly. Loading the gdal library...

> require(rgdal)

>data = readGDAL(file.choose())
/Users/andrewise/Documents/Stellies/FS 884/Shapefiles/GIS/Correct/Grabouw/30deg100.asc has GDAL driver GTiff
and has 294 rows and 250 columns
> data
Object of class SpatialGridDataFrame
Object of class SpatialGrid
Grid topology:
cellcentre.offset cellsize cells.dim
x -6561.105 100 250
y -3798028.105 100 294
SpatialPoints:
x y
[1,] -6561.10486 -3768728
[2,] -6461.10486 -3768728 ...etc

It is possible to extract the cell coordinates. However this does not produce values. The following function brings up the help file for the "raster" library.

>??raster

This library contains a number of useful analysis functions.

I assigned the variable "r" to a coarse 200x200m raster of my data.

> r <-raster(file.choose())

Checking what my variable represents:

> r

class : RasterLayer

dimensions : 147, 125, 18375 (nrow, ncol, ncell)

resolution : 200, 200 (x, y)

extent : -6611.105, 18388.9, -3798078, -3768678 (xmin, xmax, ymin, ymax)

coord. ref. : +proj=tmerc +lat_0=0 +lon_0=19 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs

values : /Users/andrewise/Documents/Stellies/FS 884/Shapefiles/GIS/Correct/Grabouw/30deg200.asc

layer name : X30deg200

The unique values of the raster:

> unique(r)
[1] 0 1 2 3 5 7 8 10 11 12 16 20 22 30 39 40 49 56 78 997 998

A list of all the raster values (many 0s, important to specify null values in the rasterization process):

> values(r)
[1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
[36] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[71] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[106] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[141] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 998 1 0 0 0 0 0 2 0 0 0 0 0 0 0 0
[176] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[211] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[246] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[281] 5 22 1 5 998 1 1 0 0 0 5 2 2 2 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

And finally, what I have been searching for, a summarized frequency table of the raster:

> freq(r)
value count
[1,] 0 15383
[2,] 1 25
[3,] 2 923
[4,] 3 1
[5,] 5 367
[6,] 7 31
[7,] 8 11
[8,] 10 159
[9,] 11 14
[10,] 12 362
[11,] 16 1
[12,] 20 80
[13,] 22 105
[14,] 30 26
[15,] 39 1
[16,] 40 22
[17,] 49 1
[18,] 56 2
[19,] 78 19
[20,] 997 38
[21,] 998 804

It seems strange that there is no plugin or similar to obtain this data. However, I would now like to somehow relate this data to the area which each value covers, to then evaluate the accuracy of the different raster resolution renderings and also with the original vector shapefile.

To visualize the matrix in R

> as.matrix(r)

[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]

[1,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA

[2,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA

[3,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA

[4,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA

[5,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA

[6,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA

Defining:

>a <- area(r)

And using the zonal function presents the following area output with area per zone which I can then use for comparison when looking at the different resolution rasters... SUCCESS!!!

> zonal(a,r,sum)
zone sum
[1,] 1 1000000
[2,] 2 36920000
[3,] 3 40000
[4,] 5 14680000
[5,] 7 1240000
[6,] 8 440000
[7,] 10 6360000
[8,] 11 560000
[9,] 12 14480000
[10,] 16 40000
[11,] 20 3200000
[12,] 22 4200000
[13,] 30 1040000
[14,] 39 40000
[15,] 40 880000
[16,] 49 40000
[17,] 56 80000
[18,] 78 760000
[19,] 997 1520000
[20,] 998 32160000

> freq(r)

value count
[1,] 1 25
[2,] 2 923
[3,] 3 1
[4,] 5 367
[5,] 7 31
[6,] 8 11
[7,] 10 159
[8,] 11 14
[9,] 12 362
[10,] 16 1
[11,] 20 80
[12,] 22 105
[13,] 30 26
[14,] 39 1
[15,] 40 22
[16,] 49 1
[17,] 56 2
[18,] 78 19
[19,] 997 38
[20,] 998 804
[21,] NA 15383

LONG WEEKEND!

# multiple layers

zonal(stack(r, r*10), z, 'sum')

Ok so I have tried multiple layers and have not been completely successful.

My steps:

Put all raster files into one folder
Identify the location of the rasters
>locat
[1] "~/Documents/Stellies/FS 884/Shapefiles/GIS/Correct/Grabouw/AllRasters"
Identify the class of each raster as "raster" individually:
>d<-raster(rastfiles[4])
> f<-raster(rastfiles[6])

> e<-raster(rastfiles[5])

> h<-raster(rastfiles[8])

> g<-raster(rastfiles[7])
Combine the frequency table of all the rasters in one matrix

> freqs<-cbind(freq(a),freq(b),freq(c),freq(d),freq(e),freq(f),freq(g),freq(h),freq(i),freq(j),freq(k),freq(l))
Convert the matrix to a data frame for easier manipulation

> freqsdf<-data.frame(freqs)
Rename the columns for easier identification

> colnames(freqs) <- c(1:24)
Create data frame again with clearer column headings

> freqsdf<-data.frame(freqs)
Extract the column names as variables

> attach(freqsdf)
Create new data frame with only one column for the types and all the other columns of the matrix values

> sumfreq<-cbind(X1,X2,X4,X6,X8,X10,X12,X14,X16,X18,X20,X22,X24)
Change the column names to desired

> colnames(sumfreqs) <- c("Type","A","B","C","D","E","F","G","H","I","J","K","L")
Export raster as comma delimited
> locatout <- "~/Documents/Stellies/FS 884/Shapefiles/GIS/Correct/Grabouw/outputfreqs.xls"

> write.table(sumfreqs,locatout, sep="\t")
Continue analysis as spreadsheet using pixel size to calculate the area.

Thursday, April 26, 2012

Raster Analysis Using R