Table of Contents
load rasters in R
all usual formats
#e.g. tiff: RasterLayer=raster("rasterLayer.tiff")
Several layers as RasterStack
All operations for RasterStacks work the same as for RasterLayers. They are performed for each layer in the stack.
HDF Files (zb MODIS)
library(MODIS) tmp=getSds("rasterLayer.hdf") RasterLayer=stack(tmp$SDS4gdal)
Raster from table or shapefile
write raster to disk
if the raster is written as rst file, it can happen that the extent changes. Then better use tif format.
Information of a raster
Generally, all important information of a raster are shown sing the print command. With the following fuctions, information can be extracted and can be writen into a variable for further usage.
extent(RasterLayer) #extent res(RasterLayer) #Aresolution proj4string(RasterLayer) #projection (sp package) crs(RasterLayer) # Projection (raster package)
# using an existing proj4 string proj4string(RasterLayer) <- CRS("+proj=utm +zone=32 +ellps=WGS84 +datum=WGS84 +units=m +no_defs") # using an EPSG code (in this case: MGI West) proj4string(RasterLayer) <- CRS("+init=epsg:31254")
#e.g with nearest neighbor method: newproj="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs" projectRaster(RasterLayer, crs=newproj,method="ngb")
Cut a rasters to a smaller extent by use of a smaller raster or an extent object
or draw the extent manually by clicking into a plot:
Show values interactively
plot(RasterLayer) click(RasterLayer) #now click on the pixels in the plot
Raster with same resolution and extent can easily calculated in the same way as you use R as a “normal” calculator:
#example: do an addition of tho rasters and multiply this with a third raster: (RasterLayer1+RasterLayer2)*RasterLayer3
#e.g use the max pixel value if both rasters overlap: mosaic(RasterLayer1,RasterLayer2,fun=max)
with rcl is a vector looking like this: c(0,10,1, 11,20,2) Here all values from 0 to including 10 will become to 1 and all values from 11 to 20 will become 2. rcl can also be a matrix with three columns: from, to, becomes.
If single values are to reclassify, you can also use a matrix of to columns: value, new value.
Using “reclassify” instead of chnaging the value is especially useful when rasterStacks are to reclassify.
Mask a Raster
All no Data value of the “mask” raster are applied on the rasterLayer
#e.g on double cell size: aggregate(RasterLayer,2)
#e.g. on half cell size: disaggregate(RasterLayer,2)
Resample to an other resolution
#use the target resolution of an other raster (Baseraster) resample(RasterLayer, Baseraster) #or to a cell size of e.g 100m by creating a base raster first Baseraster=raster(Rasterlayer,res=c(100,100)) resample(RasterLayer,Baseraster) #resample on same resolution and extent as a base raster: library(spatial.tools) spatial_sync_raster(rasterLayer, Baseraster)
For some applications, the resampling method is an important point. Per default, the method is set to 'bilinear'. For categorical raster use “ngb” nearest neighbor instead.
Unfortunately, the resampling command can take quite a long time. For large rasters you should consider to use gdal. You can use it from within R with “templateraster” is a raster which has the desired output reference system, resolution and extent (you can also specify these manually). The output is a Raster object in R but it is also written to disk as 'outputraster.tif'. Note that the inputraster is also a raster object stored on disk and not a Raster Object in R. See Usage of gdalUtils in Marburg open coursewarefor more information.
library(gdalUtils) ext<-extent(templateraster) luc<-gdalwarp('inputraster.tif', dstfile='outputraster.tif', t_srs=proj4string(templateraster), tr=res(templateraster), te=c(ext@xmin,ext@ymin,ext@xmax,ext@ymax), r="near",overwrite=TRUE,output_Raster=TRUE)
#e.g mean filter of a 5 x 5 environment focal(RasterLayer, w=matrix(1,nrow=5,ncol=5), fun=mean)
Specific filters for texture analysis from spectral data are implemented in the glcm package:
library(glcm) glcm(RasterLayer,window = c(3, 3),statistics = c("mean", "variance", "homogeneity", "contrast", "dissimilarity", "entropy", "second_moment", "correlation"))
with origin indicates the raster value from which the distance is calculated. In this case the distance is calculated from all raster cells with a value of 1. For eucliedean distance calculated from a vector object see euclidean_distance in the chapter rastervector.
Slope,Aspect, flow direction,...
#RasterLayer is a DEM terrain(RasterLayer, opt=c("slope","aspect","TPI","TRI", "roughness", "flowdir")) hillShade(terr$slope, terr$aspect,sunElevationAngle,sunAzimuthAngle)
Statistics on Patches (size, shape, etc)
This fuction comes from the gdistance package. For a cost analysis, a cost raster is required where the values contain the “cost” to cross the respective cell. Further, a start and end point is required.
costraster=RasterLayer start=SpatialPoint1 end=SpatialPoint2
Least Cost Path
Connect start and end point by the way of minimal costs:
tr=transition(costraster, function(x) 1/mean(x), directions=8) trC=geoCorrection(tr) CostPath=shortestPath(trC, start, end,output="SpatialLines")
plot(costraster) lines(CostPath) points(rbind(start,end))
For each cell, the accummulated costs to the start point are calculated.