Terrain analyses

Data preparation

Please download SRTM scenes of Sri Lanka from http://srtm.csi.cgiar.org/, unzip them and load them into R using the rgdal package.

library(rgdal)
srtm1 <- readGDAL("./data/SRTM/srtm_52_11.tif")

## ./data/SRTM/srtm_52_11.tif has GDAL driver GTiff 
## and has 6001 rows and 6001 columns

srtm2 <- readGDAL("./data/SRTM/srtm_53_11.tif")

## ./data/SRTM/srtm_53_11.tif has GDAL driver GTiff 
## and has 6001 rows and 6001 columns

str(srtm2)

## Formal class 'SpatialGridDataFrame' [package "sp"] with 4 slots
##   ..@ data       :'data.frame':  36012001 obs. of  1 variable:
##   .. ..$ band1: int [1:36012001] NA NA NA NA NA NA NA NA NA NA ...
##   ..@ grid       :Formal class 'GridTopology' [package "sp"] with 3 slots
##   .. .. ..@ cellcentre.offset: Named num [1:2] 80 5
##   .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
##   .. .. ..@ cellsize         : num [1:2] 0.000833 0.000833
##   .. .. ..@ cells.dim        : int [1:2] 6001 6001
##   ..@ bbox       : num [1:2, 1:2] 80 5 85 10
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "x" "y"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

par(mfrow = c(1,2))
image(srtm1);image(srtm2)

To conduct terrain analyses in R we use the raster package. But before we can start we need to merge the two SRTM scenes. [besides we aggregate the SRTM scene in order to accelerate computation – since these are only introductory examples]

library(raster)
srtm1 <- raster(srtm1)
srtm2 <- raster(srtm2)
srtm <- merge(srtm1,srtm2)
res(srtm)

## [1] 0.0008333333 0.0008333333

srtm <- projectRaster(srtm, res=90, crs=CRS("+init=epsg:32644"))
srtm.backup <- srtm
srtm <- aggregate(x = srtm, fact = 10, fun = mean)
srtm

## class       : RasterLayer 
## dimensions  : 622, 1235, 768170  (nrow, ncol, ncell)
## resolution  : 900, 900  (x, y)
## extent      : -166831.4, 944668.6, 552115.2, 1111915  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:32644 +proj=utm +zone=44 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /tmp/R_raster_daniel/r_tmp_2015-10-08_085734_2220_71600.grd 
## names       : layer 
## values      : -4, 2437.449  (min, max)

plot(srtm)

[please note that we could have spared two lines of code by directly importing the SRTM scenes into R using the raster package: srtm1 <- raster("./data/SRTM/srtm_52_11.tif")]

Fine, now we have a single SRTM scene - but it is very large and incorporates some unnecessary parts (India in this case and a lot of the Indian Ocean). Fortunately, we have a shapefile of your country. Let us cut the scene based on this.

library(rgdal)
sb <- readOGR(dsn = "./data/DSL250-Shp/", layer = "Boundary")

## OGR data source with driver: ESRI Shapefile 
## Source: "./data/DSL250-Shp/", layer: "Boundary"
## with 504 features
## It has 5 fields

plot(sb)

str(sb@proj4string)

## Formal class 'CRS' [package "sp"] with 1 slot
##   ..@ projargs: chr "+proj=tmerc +lat_0=7.000480277777778 +lon_0=80.77171111111112 +k=0.9999238418 +x_0=200000 +y_0=200000 +a=6377276.345 +b=6356075"| __truncated__

The plot looks as expected but be aware of differences between the SRTM and the boundary coordinate systems. Let’s see what happens when we ignore this issue

srtm <- crop(x = srtm, y = sb)
Error in .local(x, y, ...) : extents do not overlap

So, we need to reproject/transform the data first. This is fast and easy. Afterwards we can crop our SRTM scene. Last not least we store the merged, cropped, and reprojected raster as new file.

sb <- spTransform(x = sb, CRSobj = srtm@crs)
str(sb@proj4string)

## Formal class 'CRS' [package "sp"] with 1 slot
##   ..@ projargs: chr "+init=epsg:32644 +proj=utm +zone=44 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

srtm <- crop(x = srtm, y = sb)
plot(srtm)

extent(srtm)

## class       : Extent 
## xmin        : 337168.6 
## xmax        : 597268.6 
## ymin        : 654715.2 
## ymax        : 1087615

writeRaster(x = srtm, filename = "./results/srtm.tif", format = "GTiff", overwrite = TRUE)

Terrain analyses of the Kandy region using `R`

Now it is time to create some terrain parameters (you can connect R to other GIS systems to conduct all geomorphometric analyses that are available there. The only thing you need is the necessary library and the desired GIS system installed. Great options are RSAGA to conntect R to SAGA GIS and spgrass6 to connect R to GRASS GIS (works with GRASS 6 as well as GRASS 7). We will do this later. For now, we use the tools available in the raster package.

We will conduct the terrain analyses in the surroundings of Kandy on a 90x90m SRTM scene. After what we have learned so far it is straight forward to achieve such a raster.

kandy <- c(459047.6, 806090.7)
study_area <- rbind(c(kandy[1]-50000,kandy[1]+50000),
                    c(kandy[2]-50000,kandy[2]+50000))
srtm <- crop(srtm.backup, study_area)

Using the raster package we can create a lot of different terrain parameters by just one line of code. The result will be a multi-layer raster object, like you know it from multi-/hyperspectral satellite images. You can work on such scenes like you would work on ordinary vector objects in R.

srtm.tp <- terrain(x = srtm, opt = c("slope", "aspect", "TPI", "TRI", "roughness", "flowdir"), unit = "degrees", neighbors = 8)

Exercise

Plot the TPI raster

What kind of object did you create using the terrain function above?

What is the cell size of the object?

What is the mean slope for the raster?

What does it mean to use 8 neighbors?

You should be familiar with slope and aspect but do you know the what the other parameter are? A look in the help offers insights (they are based on (???)):

“TRI (Terrain Ruggedness Index) is the mean of the absolute differences between the value of a cell and the value of its 8 surrounding cells.
TPI (Topographic Position Index) is the difference between the value of a cell and the mean value of its 8 surrounding cells.
Roughness is the difference between the maximum and the minimum value of a cell and its 8 surrounding cells" (from ?raster::terrain)

The help also explains that we can use focal functions in order the adapt the approaches to our needs. A focal function corresponds to a moving window. It uses a matrix of weights for the neighborhood of the focal cells.

TPI for different neighborhood size:

tpiw <- function(x, w=5) {
    m <- matrix(1/(w^2-1), nc=w, nr=w)
    m[ceiling(0.5 * length(m))] <- 0 # set the centre cell 0
    f <- focal(x, m, pad = TRUE) # apply moving window
    x - f
}
tpi15 <- tpiw(x = srtm, w=15)
tpi31 <- tpiw(x = srtm, w=31)
#par(mfrow=c(1,3))
#plot(srtm.tp$tpi)
#plot(tpi15)
#plot(tpi31)

For better visualisation we calculate a hillshade and overlay the TPI raster.

srtm.hs <- hillShade(slope = terrain(srtm, opt="slope"),
    aspect = terrain(srtm, opt="aspect"),
    angle= 150, direction = 45, normalize = TRUE)
plot(srtm.hs,
    col = gray.colors(n= 255, start=.2, end=.9),
    legend = FALSE)
plot(tpi31,
    col=colorRampPalette(c("red", "white", "blue"))(255),
    alpha = .5, add=TRUE)

Excercise

TPI is scale dependent. What does it mean? How can the following plot be used to explain the phenomenon?

library(rasterVis)
#srtmTheme <- rasterTheme(region=terrain.colors(200))
#levelplot(srtm, par.settings = srtmTheme)
levelplot(brick(tpi31, tpi15,srtm.tp$tpi), par.settings = RdBuTheme)

Create elevation profile

a <- Line(coords = cbind(c(kandy[1]-50000,kandy[1]+50000),c(kandy[2] ,kandy[2])))
b <- Lines(list(a),ID="1")
c <- SpatialLines(list(b), proj4string=srtm@crs)
e.prof <- extract(x = srtm, y = c)
t.prof <- extract(x = srtm.tp, y = c)

prof.x <- cellFromLine(srtm, c)
prof.x <- xyFromCell(srtm,prof.x[[1]])

plot(srtm)
lines(c)

plot(prof.x[,1],e.prof[[1]], type = "l",lwd=.6, ylab="Elevation (m)", xlab="UTM coordinate in Meter")
abline(v = kandy, col = "red")

#plot(t.prof[[1]][,2], type = "l",lwd=.2, ylab="TPI")

Terrain parameter of the Kandy region using GRASS GIS

Be aware, for the following to work you need to have GRASS GIS installed on your system.

Load library, initialize the GRASS environment with default values (everything is stored as temporary files; you can change this if you want to come back to earlier stages of your work)

library(spgrass6)
loc <-  initGRASS("/usr/local/grass-7.0.0svn/", home=tempdir() ,mapset = "PERMANENT", override = TRUE)

Now we check and set our spatial reference and load our SRTM from R into GRASS.

## check the location and define it according to our data
execGRASS("g.proj", flags = c("p"))

## XY location (unprojected)

execGRASS("g.proj", flags = c("c"), parameters = list(proj4="+init=epsg:32644"))

## Default region was updated to the new projection, but if you have multiple
## mapsets `g.region -d` should be run in each to update the region from the
## default
## Projection information updated

## load data from R to grass
library(rgdal)
writeRAST6(x = as(srtm, "SpatialGridDataFrame"), vname="dem", flags = c("overwrite"))

Let us inspect the dataset and the characteristics of our GRASS region.

## GRASS - show raster in mapset
execGRASS("r.info", parameters = list(map = "dem"))

##  +----------------------------------------------------------------------------+
##  | Map:      dem                            Date: Thu Oct  8 08:57:49 2015    |
##  | Mapset:   PERMANENT                      Login of Creator: daniel          |
##  | Location: file8ac18ffbab0                                                  |
##  | DataBase: /tmp/RtmpN7rsL4                                                  |
##  | Title:     ( dem )                                                         |
##  | Timestamp: none                                                            |
##  |----------------------------------------------------------------------------|
##  |                                                                            |
##  |   Type of Map:  raster               Number of Categories: 0               |
##  |   Data Type:    FCELL                                                      |
##  |   Rows:         1112                                                       |
##  |   Columns:      1111                                                       |
##  |   Total Cells:  1235432                                                    |
##  |        Projection: UTM (zone 44)                                           |
##  |            N: 856135.18883952    S: 756055.18883952   Res:    90           |
##  |            E: 509068.57873155    W: 409078.57873155   Res:    90           |
##  |   Range of data:    min = 16.1446  max = 2514.008                          |
##  |                                                                            |
##  |   Data Description:                                                        |
##  |    generated by r.in.gdal                                                  |
##  |                                                                            |
##  |   Comments:                                                                |
##  |    r.in.gdal input="/tmp/RtmpN7rsL4/file8ac18ffbab0/PERMANENT/.tmp/Uni/\   |
##  |    X839" output="dem" memory=300 offset=0                                  |
##  |                                                                            |
##  +----------------------------------------------------------------------------+

execGRASS("g.region", flags = c("p"))

## projection: 1 (UTM)
## zone:       44
## datum:      wgs84
## ellipsoid:  wgs84
## north:      1
## south:      0
## west:       0
## east:       1
## nsres:      1
## ewres:      1
## rows:       1
## cols:       1
## cells:      1

As you can see there is a difference in the extent and resolution of our GRASS area and the SRTM. Let us fix this.

## adjust Regions resolution!
execGRASS("g.region", parameters = list(raster = "dem",res = as.character(res(srtm)[1])))
execGRASS("g.region", flags = c("p"))

## projection: 1 (UTM)
## zone:       44
## datum:      wgs84
## ellipsoid:  wgs84
## north:      856135.18883952
## south:      756055.18883952
## west:       409078.57873155
## east:       509068.57873155
## nsres:      90
## ewres:      90
## rows:       1112
## cols:       1111
## cells:      1235432

After fixing this we will now create differnt terrain paramters.

Example 1: Curvature

parseGRASS("r.param.scale")

execGRASS("r.param.scale", parameters = list(input = "dem",
                               output = "profcurv",
                               method = "profc",
                               size = 13),
          flags = c("overwrite"))

profcurv <- raster(readRAST6("profcurv"))
plot(srtm.hs,
    col = gray.colors(n= 255, start=.2, end=.9),
    legend = FALSE)
plot(profcurv, alpha=.5, add=TRUE)
title("Profile Curvature")

Example 2:Simple landform classification

execGRASS("r.param.scale", parameters = list(input = "dem",
                               output = "landforms",
                               method = "feature",
                               size = 15),
          flags = c("overwrite"))

landforms <- raster(readRAST6("landforms"))

landforms <- ratify(landforms)
tmp <- levels(landforms)[[1]]
tmp$Landform <- c("Planar","Pit","Channel","Pass","Ridge","Peak")
levels(landforms) <- tmp
levels(landforms)

#library(rasterVis)
#levelplot(landforms)

plot(srtm.hs,
    col = gray.colors(n= 255, start=.2, end=.9),
    legend = FALSE)
plot(landforms, col=colorRampPalette(c("white", "darkblue","blue", "yellow","brown","black"))(6),alpha=.5, add=TRUE, legend = FALSE)
title("Landforms")
mycol <- colorRampPalette(c("white", "darkblue","blue", "yellow","brown","black"))(6)
legend(x="topright",legend = tmp$Landform, fill = mycol)
points(SpatialPoints(coords=cbind(kandy[1],kandy[2]), proj4string = srtm@crs), pch=19)

Exercise

What does the different terrain parameters show? How are they characterized? Answer the question using the help of the corresponding GRASS GIS package.

Watershed delineation

The GRASS packages r.watershed is a great tool for fast and comprehensive analyses of basic terrain-hydrological analyses. We will create different rasters that can be used for subsequent and more adavanced analyses.

Besides r.watershed there is the r.stream.* family that has complementary and besides this many additional features, e.g. automatic extraction of stream orders and calculation of stream and catchment characteristics. You need to install the package in GRASS first using g.extension. For more information have a look at https://grasswiki.osgeo.org/wiki/R.stream.*_modules

## catchment analyses
parseGRASS("r.watershed")

execGRASS("r.watershed", parameters = list(elevation = "dem",
                             threshold = 15000,
                             accumulation = "fac",
                             basin = "basins",
                             drainage = "fdir",
                             stream = "streams_r"
                                           ),
          flags = c("overwrite","b","a"))

## Warning in execGRASS("r.watershed", parameters = list(elevation = "dem", : The command:
## r.watershed --overwrite -b -a elevation=dem threshold=15000 accumulation=fac basin=basins drainage=fdir stream=streams_r
## produced at least one warning during execution:
## SECTION 1a (of 5): Initiating Memory.
## SECTION 1b (of 5): Determining Offmap Flow.
##    0%   4%   8%  12%  16%  20%  24%  28%  32%  36%  40%  44%  48%  52%  56%  60%  64%  68%  72%  76%  80%  84%  88%  92%  96% 100%
## SECTION 2: A* Search.
##    0%   2%   4%   6%   8%  10%  12%  14%  16%  18%  20%  22%  24%  26%  28%  30%  32%  34%  36%  38%  40%  42%  44%  46%  48%  50%  52%  54%  56%  58%  60%  62%  64%  66%  68%  70%  72%  74%  76%  78%  80%  82%  84%  86%  88%  90%  92%  94%  96%  98% 100%
## SECTION 3a: Accumulating Surface Flow with MFD.
##    1%   3%   5%   7%   9%  11%  13%  15%  17%  19%  21%  23%  25%  27%  29%  31%  33%  35%  37%  39%  41%  43%  45%  47%  49%  51%  53%  55%  57%  59%  61%  63%  65%  67%  69%  71%  73%  75%  77%  79%  81%  83%  85%  87%  89%  91%  93%  95%  97%  99% 100%
## SECTION 3b: Adjusting drainage directions.
##    1%   3%   5%   7%   9%  11%  13%  15%  17%  19%  21%  23%  25%  27%  29%  31%  33%  35%  37%  39%  41%  43%  45%  47%  49%  51%  53%  55%  57%  59%  61%  63%  65%  67%  69%  71%  73%  75%  77%  79%  81%  83%  85%  87%  89%  91%  93%  95%  97%  99% 100%
## SECTION 4: Watershed determination.
##    0%   2%   4%   6%   8%  10%  12%  14%  16%  18%  20%  22%  24%  26%  28%  30%  32%  34%  36%  38%  40%  42%  44%  46%  48%  50%  52%  54%  56%  58%  60%  62%  64%  66%  68%  70%  72%  74%  76%  78%  80%  82%  84%  86%  88%  90%  92%  94%  96%  98% 100%
## SECTION 5: Closing Maps.
##  100%
## WARNING: Writing out only positive flow accumulation values.
## WARNING: Cells with a likely underestimate for flow accumulation can no
##          longer be identified.

fac <- raster(readRAST6("fac"))
plot(fac)

fac2 <- stretch(fac, minq=.1,maxq=.9)

plot(srtm.hs,
    col = gray.colors(n= 255, start=.2, end=.9),
    legend = FALSE)
plot(fac2, col = colorRampPalette(c("yellow","green","blue"))(255), alpha = .4, add=TRUE)
title("Flow accumulation - 8bit scaled")

fdir <- raster(readRAST6("fdir"))
#plot(fdir, col = gray.colors(255,0,.6))

basins <- raster(readRAST6("basins"))
plot(basins, col = rainbow(50), legend = FALSE)

parseGRASS("r.to.vect")

execGRASS("r.to.vect", parameters = list(input = "streams_r",
                           output = "streams",
                           type = "line"
                           ),
          flags = c("v","overwrite"))          

execGRASS("v.out.ogr", parameters = list(input = "streams",
                           type = "line",
                           output = "/media/daniel/homebay/Teachings/WS_2015-2016/Statistics_SriLanka/results/",
                           format = "ESRI_Shapefile"
                           ),
          flags = c("overwrite"))

## Warning in execGRASS("v.out.ogr", parameters = list(input = "streams", type = "line", : The command:
## v.out.ogr --overwrite input=streams type=line output=/media/daniel/homebay/Teachings/WS_2015-2016/Statistics_SriLanka/results/ format=ESRI_Shapefile
## produced at least one warning during execution:
## WARNING: OGR layer <streams> already exists and will be overwritten
## Exporting 157 features...
##    5%  11%  17%  23%  29%  35%  41%  47%  53%  59%  65%  71%  77%  83%  89%  95% 100%
## v.out.ogr complete. 157 features (Line String type) written to <streams>
## (ESRI_Shapefile format).

streams <- readOGR("./results","streams")

plot(srtm.hs,
    col = gray.colors(n= 255, start=.2, end=.9),
    legend = FALSE)
plot(basins, col = rainbow(30), legend = FALSE, alpha = .3, add=TRUE)
lines(streams, col="blue")
points(SpatialPoints(coords=cbind(kandy[1],kandy[2]), proj4string = srtm@crs), pch=19)

## optional: add the "real" Mahaweli Ganga
# tmp<-readOGR("./data/DSL250-Shp", "Streams")
# tmp <- spTransform(x = tmp, CRSobj = srtm@crs)
# plot(tmp[tmp@data$NAME=="Mahaweli Ganga",], add=TRUE)

title("Basins and streams")

Landform classification

Landform classification based on machine learning ideas and principal component analysis.

First we will add some terrain parameters calculated with GRASS to our raster brick object srtm.tp.

srtm.tp$profcurv <- brick(profcurv)
srtm.tp$tpi31 <- brick(tpi31)
srtm.tp

## class       : RasterStack 
## dimensions  : 1112, 1111, 1235432, 8  (nrow, ncol, ncell, nlayers)
## resolution  : 90, 90  (x, y)
## extent      : 409078.6, 509068.6, 756055.2, 856135.2  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:32644 +proj=utm +zone=44 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## names       :           tri,           tpi,     roughness,         slope,        aspect,       flowdir,      profcurv,         tpi31 
## min values  :  0.000000e+00, -8.923577e+01,  0.000000e+00,  0.000000e+00,  0.000000e+00,  1.000000e+00, -1.552875e-03, -3.479182e+02 
## max values  :  1.709198e+02,  8.886958e+01,  4.968500e+02,  6.712930e+01,  3.599997e+02,  1.280000e+02,  1.965019e-03,  5.139620e+02

tp <- as(srtm.tp, "SpatialGridDataFrame")

# principal component analysis
pc.dem <- prcomp(~tri+tpi+roughness+slope+flowdir+profcurv+tpi31, data = tp@data, scale=TRUE,na.action = na.exclude) # na.exclude and not na.omit, since it is intended to keep NAs (because these represent the ocean and of course lacks in the data)
biplot(pc.dem, arrow.len = 0.1, xlabs = rep(".", length(pc.dem$x[,1])),main = "PCA biplot") # takes time...

# in case the biplot shows that some variables are correlated
#pc.dem <- prcomp(~TRI+TWI+TPI, data = dem.param@data, scale=TRUE,na.action = na.exclude) # na.exclude and not na.omit, since it is intended to keep NAs (because these represent the ocean and of course lacks in the data)
#biplot(pc.dem, arrow.len = 0.1, xlabs = rep(".", length(pc.dem$x[,1])),main = "PCA biplot") # takes time...

demdata <- as.data.frame(pc.dem$x)

tmp<-na.omit(demdata)

wss <- (nrow(tmp)-1)*sum(apply(tmp,2,var))
for (i in 2:20) {
     wss[i] <- sum(kmeans(tmp,centers=i)$withinss)
}

## Warning: Quick-TRANSfer stage steps exceeded maximum (= 58482100)

## Warning: Quick-TRANSfer stage steps exceeded maximum (= 58482100)

## Warning: Quick-TRANSfer stage steps exceeded maximum (= 58482100)

## Warning: Quick-TRANSfer stage steps exceeded maximum (= 58482100)

## Warning: Quick-TRANSfer stage steps exceeded maximum (= 58482100)

## Warning: Quick-TRANSfer stage steps exceeded maximum (= 58482100)

## Warning: Quick-TRANSfer stage steps exceeded maximum (= 58482100)

## Warning: Quick-TRANSfer stage steps exceeded maximum (= 58482100)

## Warning: Quick-TRANSfer stage steps exceeded maximum (= 58482100)

## Warning: Quick-TRANSfer stage steps exceeded maximum (= 58482100)

## Warning: Quick-TRANSfer stage steps exceeded maximum (= 58482100)

## Warning: Quick-TRANSfer stage steps exceeded maximum (= 58482100)

plot(wss, type = "l")

kmeans.dem <- kmeans(na.omit(demdata),centers=5,iter.max=500) #kmeans 
# two small functions to implement the omitted NAs in the end
demdata$cluster <- rep(NA, nrow(demdata))
demdata[names(kmeans.dem$cluster), "cluster"] <- kmeans.dem$cluster
tp$kmeans.dem <- demdata$cluster
#dem.param$landform <- as.factor(demdata$cluster)
#summary(dem$landform)
#tp@data$kmeans.dem <- as.numeric(dem.param@data$kmeans.dem) # for export, gdal needs numeric values

# remove all bands except the fuzzy kmeans results
lc <- tp
lc@data <- subset(lc@data, select = kmeans.dem) # just select the relevant data
## another possibility:
## lc <- raster(lc,5) # create raster object using band 5

plot(raster(lc))

Geomorphometric Analyses

ESR 512 - GIS in Geostatistics

Postgraduate Institute of Science

University of Peradeniya

Thusitha Wagalawatta & Daniel Knitter

10/2015

Terrain analyses

Data preparation

Terrain analyses of the Kandy region using `R`

TPI for different neighborhood size:

Terrain parameter of the Kandy region using GRASS GIS

Example 1: Curvature

Example 2:Simple landform classification

Watershed delineation

Landform classification

Geomorphometric Analyses

ESR 512 - GIS in Geostatistics

Postgraduate Institute of Science

University of Peradeniya

Thusitha Wagalawatta & Daniel Knitter

10/2015

Terrain analyses

Data preparation

Terrain analyses of the Kandy region using R

TPI for different neighborhood size:

Terrain parameter of the Kandy region using GRASS GIS

Example 1: Curvature

Example 2:Simple landform classification

Watershed delineation

Landform classification

Terrain analyses of the Kandy region using `R`