User Tools

Site Tools


r:r-gis:raster

1-Raster Operations

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.

  RasterStack=stack("rasterLayer1.tiff",
    "rasterLayer2.tiff","rasterLayer3.tiff")

HDF Files (zb MODIS)

  library(MODIS)
  tmp=getSds("rasterLayer.hdf")
  RasterLayer=stack(tmp$SDS4gdal)

Raster from table or shapefile

write raster to disk

  writeRaster(rasterLayer, filename="filename.tif")

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)

Assign/define projection

  # 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")  

Project raster

  #e.g with nearest neighbor method:
  newproj="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
  projectRaster(RasterLayer, crs=newproj,method="ngb")

Crop

Cut a rasters to a smaller extent by use of a smaller raster or an extent object

  crop(raster1,smallerRaster/extentObject)

or draw the extent manually by clicking into a plot:

  crop(raster1,drawExtent())

Show values interactively

  plot(RasterLayer)
  click(RasterLayer)
  #now click on the pixels in the plot

Raster calculator

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

Mosaic

  #e.g use the max pixel value if both rasters overlap:
  mosaic(RasterLayer1,RasterLayer2,fun=max)

Reclassify

  rasterLayer[rasterLayer>/</==/!=oldValue]=newValue

or

reclassify(rasterLayer, rcl)

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

  mask(rasterLayer, mask)

Aggregation

  #e.g on double cell size:
  aggregate(RasterLayer,2)

Disaggregation

  #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) 

Filter

#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"))

Euclidean distance

gridDistance(RasterLayer,origin=1)

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 3-Raster-Vector Operations.

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)

Identify Patches

patches<-clump(rasterLayer)

Statistics on Patches (size, shape, etc)

library(SDMTools)
PatchStat(patches)

Cost analysis

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")

Show results:

  plot(costraster)
  lines(CostPath)
  points(rbind(start,end))

Akkumulate costs

For each cell, the accummulated costs to the start point are calculated.

  accCost(trC,start)

Hanna Meyer

r/r-gis/raster.txt · Last modified: 2018/12/23 19:46 (external edit)