-
Notifications
You must be signed in to change notification settings - Fork 0
/
ktpi.R
executable file
·275 lines (239 loc) · 19.2 KB
/
ktpi.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
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
#!/usr/bin/env Rscript
suppressMessages(library(docopt))
suppressMessages(library(jsonlite))
# ktpi.R neighbours -c 3 -r 4 -w 2 -x 4 -y 3 -z 5 -m 1 -n 500 -k 100 -l tiles.txt
# ktpi.R neighbours -p features/ -c 3 -r 4 -a .tif -w 2 -x 4 -y 3 -z 5 -m 1 -n 500 -k 100 -l tiles.txt
# ktpi.R statistic -p features/ -q dems/ -c 3 -r 4 -a .tif -w 1 -x 99 -y 1 -z 99 -m 1 -n 500 -d 5 -l tiles.txt -u output/
# ktpi.R terrain -p features/ -q dems/ -c 3 -r 4 -a .tif -w 1 -x 99 -y 1 -z 99 -m 1 -n 500 -d 5 -l tiles.txt -u output/
# ktpi.R ktpi -p features/ -q dems/ -c 3 -r 4 -a .tif -w 1 -x 99 -y 1 -z 99 -m 1 -n 500 -d 5 -k 100 -l tiles.txt -u output/
# ktpi.R kaspSlp -p features/ -q dems/ -c 3 -r 4 -a .tif -w 1 -x 99 -y 1 -z 99 -m 1 -n 500 -d 5 -k 100 -o across -l tiles.txt -u output/
# ktpi.R clicommands --ktpi-function statistic --ktpi-function terrain -p features/ -q dems/ -a .tif -w 1 -x 99 -y 1 -z 99 -m 1 -n 500 -d 1 -d 5 -d 10 -d 20 -l tiles.txt -u output/
# ktpi.R clicommands --ktpi-function ktpi -p features/ -q dems/ -a .tif -w 1 -x 99 -y 1 -z 99 -m 1 -n 500 -d 5 -f 20 -t 50 -s 10 -d 10 -f 50 -t 200 -s 50 -l tiles.txt -u output/
# ktpi.R clicommands --ktpi-function kaspSlp -p features/ -q dems/ -a .tif -w 1 -x 99 -y 1 -z 99 -m 1 -n 500 -d 5 -f 20 -t 50 -s 10 -d 10 -f 50 -t 200 -s 50 -o uphill -o across -o downhill -l tiles.txt -u output/
# ktpi.R jsonmessages --ktpi-function statistic --ktpi-function terrain -p features/ -q dems/ -a .tif -w 1 -x 99 -y 1 -z 99 -m 1 -n 500 -d 1 -d 5 -d 10 -d 20 -l tiles.txt -u output/
# ktpi.R jsonmessages --ktpi-function ktpi -p features/ -q dems/ -a .tif -w 1 -x 99 -y 1 -z 99 -m 1 -n 500 -d 5 -f 20 -t 50 -s 10 -d 10 -f 50 -t 200 -s 50 -l tiles.txt -u output/
# ktpi.R jsonmessages --ktpi-function kaspSlp -p features/ -q dems/ -a .tif -w 1 -x 99 -y 1 -z 99 -m 1 -n 500 -d 5 -f 20 -t 50 -s 10 -d 10 -f 50 -t 200 -s 50 -o uphill -o across -o downhill -l tiles.txt -u output/
# ktpi.R jsonmessages --ktpi-function statistic --ktpi-function terrain --ktpi-function ktpi --ktpi-function kaspSlp --ktpi-function kaspDir --ktpi-function kaspSDir --ktpi-function kaspCDir --ktpi-function kaspSlpSDir --ktpi-function kaspSlpCDir --ktpi-function kaspSlpEle2 --ktpi-function kaspSlpEle2SDir --ktpi-function kaspSlpEle2CDir --ktpi-function kaspSlpLnEle --ktpi-function kaspSlpLnEleSlpSDir --ktpi-function kaspSlpLnEleSlpCDir -p features/ -q dems/ -a .tif -w 1 -x 3 -y 1 -z 2 -m 1 -n 4000 -d 1 -f 10 -t 20 -s 5 -d 5 -f 20 -t 100 -s 20 -d 10 -f 100 -t 200 -s 20 -d 10 -f 250 -t 500 -s 50 -d 10 -f 600 -t 1000 -s 100 -d 20 -f 1000 -t 2000 -s 200 -d 20 -f 2250 -t 2500 -s 250 -o across -o uphill -o downhill -u output/ > json.txt
# ktpi.R jsonmessages --ktpi-function statistic --ktpi-function terrain --ktpi-function ktpi --ktpi-function kaspSlp --ktpi-function kaspDir --ktpi-function kaspSDir --ktpi-function kaspCDir --ktpi-function kaspSlpSDir --ktpi-function kaspSlpCDir --ktpi-function kaspSlpEle2 --ktpi-function kaspSlpEle2SDir --ktpi-function kaspSlpEle2CDir --ktpi-function kaspSlpLnEle --ktpi-function kaspSlpLnEleSlpSDir --ktpi-function kaspSlpLnEleSlpCDir -p features/ -q dems/ -a .tif -w 1 -x 14 -y 1 -z 11 -m 1 -n 2000 -d 1 -f 5 -t 20 -s 5 -d 5 -f 20 -t 60 -s 10 -d 10 -f 60 -t 200 -s 20 -d 20 -f 200 -t 600 -s 50 -d 20 -f 600 -t 1000 -s 100 -d 40 -f 1000 -t 2000 -s 200 -o across -o uphill -o downhill -l limit-tiles.txt -u output/ > kc_stat_terr_ktpi_kasp.json
'Usage:
ktpi.R (statistic | terrain) (--folder <fpre>) (--dem-folder <dpre>) (--tile-col <tcol>) (--tile-row <trow>) (--extension <ext>) (--tile-col-min <cmin>) (--tile-col-max <cmax>) (--tile-row-min <rmin>) (--tile-row-max <rmax>) (--raster-cell-size <csize>) (--raster-cells <rcell>) [--dem-calc-size <dsize>] [--limit-tiles <tiles-csv>] (--output-folder <ofolder>) [--exp-rast]
ktpi.R ktpi (--folder <fpre>) (--dem-folder <dpre>) (--tile-col <tcol>) (--tile-row <trow>) (--extension <ext>) (--tile-col-min <cmin>) (--tile-col-max <cmax>) (--tile-row-min <rmin>) (--tile-row-max <rmax>) (--raster-cell-size <csize>) (--raster-cells <rcell>) [--dem-calc-size <dsize>] [--kernel-size <ksize>] [--limit-tiles <tiles-csv>] (--output-folder <ofolder>) [--exp-rast]
ktpi.R (kaspSlp | kaspDir | kaspSDir | kaspCDir | kaspSlpSDir | kaspSlpCDir | kaspSlpEle2 | kaspSlpEle2SDir | kaspSlpEle2CDir | kaspSlpLnEle | kaspSlpLnEleSlpSDir | kaspSlpLnEleSlpCDir) (--folder <fpre>) (--dem-folder <dpre>) (--tile-col <tcol>) (--tile-row <trow>) (--extension <ext>) (--tile-col-min <cmin>) (--tile-col-max <cmax>) (--tile-row-min <rmin>) (--tile-row-max <rmax>) (--raster-cell-size <csize>) (--raster-cells <rcell>) [--dem-calc-size <dsize>] [--kernel-size <ksize>] [--orientation <orient>] (--output-folder <ofolder>) [-l|--limit-tiles <tiles-csv>] [--exp-rast]
ktpi.R neighbours [--folder <pre>] (-c <fcol>) (-r <frow>) [--extension <app>] (--tile-col-min <cmin>) (--tile-col-max <cmax>) (--tile-row-min <rmin>) (--tile-row-max <rmax>) (--raster-cells <rcell>) (--raster-cell-size <csize>) (-k <ksize>) [-l|--limit-tiles <tiles-csv>]
ktpi.R clicommands (--ktpi-function <ktpi-func>)... (--folder <fpre>) (--dem-folder <dpre>) (--extension <app>) (--tile-col-min <cmin>) (--tile-col-max <cmax>) (--tile-row-min <rmin>) (--tile-row-max <rmax>) (--raster-cells <rcell>) (--raster-cell-size <csize>) (-d <dsize>... [-f <kfrom> -t <kto> -s <kstep>]...) [-o|--orientation <orient>...] [-l|--limit-tiles <tiles-csv>] (--output-folder <ofolder>) [--exp-rast]
ktpi.R jsonmessages (--ktpi-function <ktpi-func>)... (--folder <fpre>) (--dem-folder <dpre>) (--extension <app>) (--tile-col-min <cmin>) (--tile-col-max <cmax>) (--tile-row-min <rmin>) (--tile-row-max <rmax>) (--raster-cells <rcell>) (--raster-cell-size <csize>) (-d <dsize>... [-f <kfrom> -t <kto> -s <kstep>]...) [-o|--orientation <orient>...] [-l|--limit-tiles <tiles-csv>] (--output-folder <ofolder>) [--exp-rast]
ktpi.R [-h|--help]
kpti.R --version
ktpi.R info
Options:
-h --help show this screen.
--version show version.
statistic runs statistic indices at the defined cell size. ie. statistic dsize=5, will develop statistic indices at the 5m cell size.
terrain runs terrain indices at the defined cell size. ie. terrain dsize=10, will develop terrain indices at the 10m cell size.
ktpi runs topographic position indices at the defined cell and kernel neighbourhood size. ie. ktpi dsize=20 kernel-size=800, will develop topographic position indices (sd & mean_diff) for a 20m dem cell size at the 800m kernel neighbourhood.
kaspSlp runs slope indices at the defined orientation (across, uphill, downhill) at the defined cell and kernel neighbourhood size. ie. ktpi-aspect ktpislp uphill dsize=10 kernel-size=500, will develop slope indices for a 10m dem cell size at the 500m kernel ring for DEM values uphill of the kernel centroid.
kaspDir runs direction indices at the defined orientation (across, uphill, downhill) at the defined cell and kernel neighbourhood size.
kaspSDir runs sin direction indices at the defined orientation (across, uphill, downhill) at the defined cell and kernel neighbourhood size.
kaspCDir runs cos direction indices at the defined orientation (across, uphill, downhill) at the defined cell and kernel neighbourhood size.
kaspSlpSDir runs Al Stage slope * sin direction
kaspSlpCDir runs Al Stage slope * cos direction
kaspSlpEle2 runs Al Stage slope * dem^2
kaspSlpEle2SDir runs Al Stage slope * dem^2 * sin direction
kaspSlpEle2CDir runs Al Stage slope * dem^2 * sin direction
kaspSlpLnEle runs Al Stage slope * natural log of dem
kaspSlpLnEleSlpSDir runs Al Stage slope * natural log of dem * slp * sin direction
kaspSlpLnEleSlpCDir runs Al Stage slope * natural log of dem * slp * cos direction
neighbours gets a list of the neighbouring tiles of a given raster tile, based on the given raster dimensions, and a given kernel size, and a text file of tiles to filter (column of "col/row").
clicommands creates a list of CLI commands on min/max column/row, dem calculation sizes, kernel from-to-step, orientation, and a text file of tiles to filter
jsonmessages creates a list of JSON messages on min/max column/row, dem calculation sizes, kernel from-to-step, orientation, and a text file of tiles to filter
-p <fpre>, --folder <fpre> prepended input feature raster source folder (/<folder>/)
-q <dpre>, --dem-folder <dpre> prepended input DEM raster source folder (/<folder>/)
-c <fcol>, --tile-col <fcol> tile column
-r <frow>, --tile-row <frow> tile row
-a <ext>, --extension <ext> appended tile file extension (.tif)
-w <cmin>, --tile-col-min <cmin> project tile min column value
-x <cmax>, --tile-col-max <cmax> project tile max column value
-y <rmin>, --tile-row-min <rmin> project tile min row value
-z <rmax>, --tile-row-max <rmax> project tile max row value
-m <csize>, --raster-cell-size <csize> raster cell size (m)
-n <rcell>, --raster-cells <rcell> raster tile size (cells)
-d <dsize>, --dem-calc-size <dsize> DEM cell size to calculate indices. ie. 5 = 5m [DEFAULT = raster cell size, maximum = cell size * tile cells] ( recalculated to be nearest evenly divisible into tile size: for (i in 1:(tilecells)) {if(tilecells%%i == 0) {print(i)}} ).
-k <ksize>, --kernel-size <ksize> kernel neighbourhood size, in ground units, to calculate ktpi indices. ie. 50 = 50m (DEFAULT = demCalcSize x 5)
-f <kfrom>, --kernel-from <kfrom> kernel from range
-t <kto>, --kernel-to <kto> kernel to ranged
-s <kstep>, --kernel-step <kstep> kernel range step
-o <orient>, --orientation <orient> ktpi-aspect kernel orientation: across, uphill, downhill
-l <tiles>, --limit-tiles <tiles> text file of tiles to filter (column of "col/row") [default: none]
--ktpi-function <ktpi-func> ktpi functions to generate CLI/JSON arguments
-u <ofolder>, --output-folder <ofolder> existing output data folder (/<folder>/)
--exp-rast exports indices rasters.
' -> doc
args <- docopt(doc)
inputDataPath <- file.path(Sys.getenv('HRIS_DATA'), 'input')
outputDataPath <- file.path(Sys.getenv('HRIS_DATA'), 'output')
dir.create(outputDataPath, showWarnings = FALSE)
wd <- Sys.getenv('HRIS_DATA')
setwd(wd)
options(warn = -1)
source_local <- function(fname){
argv <- commandArgs(trailingOnly = FALSE)
base_dir <- dirname(substring(argv[grep("--file=", argv)], 8))
source(paste(base_dir, fname, sep="/"))
}
ensureDataPath <- function(dataPath){
if (!startsWith(dataPath, "input")) {
dataPath <- file.path("input", dataPath)
}
return(dataPath)
}
if (args$'info') {
print(.libPaths())
print(sessionInfo())
print(version)
}
tilesfile <- args$'limit-tiles'
tiles <- "none"
if (tilesfile != "none") {
tilesFileFullPath <- file.path(Sys.getenv('HRIS_DATA'), ensureDataPath(tilesfile))
tiles <- read.csv(tilesFileFullPath, stringsAsFactors=FALSE, header=FALSE)$V1
}
if (args$'statistic' | args$'terrain' | args$'ktpi' | args$'kaspSlp' | args$'kaspDir' |
args$'kaspSDir' | args$'kaspCDir' | args$'kaspSlpSDir' | args$'kaspSlpCDir' |
args$'kaspSlpEle2' | args$'kaspSlpEle2SDir' | args$'kaspSlpEle2CDir' |
args$'kaspSlpLnEle' | args$'kaspSlpLnEleSlpSDir' | args$'kaspSlpLnEleSlpCDir') {
source_local("lib/ktpi_util.R")
source_local("lib/ktpi_aspect.R")
source_local("lib/ktpi_terrain.R")
folderPath <- ensureDataPath(args$'folder')
demFolderPath <- ensureDataPath(args$'dem-folder')
# test the actual image dimension matches the given image dimension (-n / --raster-cells)
imagePath <- paste(folderPath, args$'tile-col', "/", args$'tile-row', args$'extension', sep = "")
actualRasterCells <- ncol(raster(imagePath))
if (actualRasterCells != args$'raster-cells') {
stop(paste("Your image size is ", actualRasterCells, " please change your raster-cells (-n) value.", sep = ""))
}
# gets feature neighbouring tile files based on the kernel neighbourhood size
neighbourFileList <- getNeighbours(folderPath, args$'tile-col', args$'tile-row', args$'extension',
args$'tile-col-min', args$'tile-col-max', args$'tile-row-min', args$'tile-row-max',
args$'raster-cells', args$'raster-cell-size', args$'kernel-size', tiles)
print(neighbourFileList)
# merges feature neighbouring tile files
print ("merges feature neighbouring tile files")
featureNeighbourRaster <- mergeRasters(neighbourFileList)
# gets dem neighbouring tile files based on the kernel neighbourhood size
neighbourFileList <- getNeighbours(demFolderPath, args$'tile-col', args$'tile-row', args$'extension',
args$'tile-col-min', args$'tile-col-max', args$'tile-row-min', args$'tile-row-max',
args$'raster-cells', args$'raster-cell-size', args$'kernel-size', "none")
print(neighbourFileList)
# merges dem neighbouring tile files
print ("merges dem neighbouring tile files")
demNeighbourRaster <- mergeRasters(neighbourFileList)
# ensure featureNeighbourRaster & demNeighbourRasterExtent are equivalent
featureNeighbourRaster <- extend(featureNeighbourRaster, extent(demNeighbourRaster), value=NA)
# gets initial indice table with unique trFeatId and cell counts
indic <- getFeatureIdCount(folderPath,
args$'tile-col',
args$'tile-row',
args$'extension')
indicFeatCountInit <- nrow(indic)
# runs statistic indices
if (args$'statistic') {
ktpiFunction <- "statistic"
featStat <- statisticIndices(folderPath, args$'tile-col', args$'tile-row', args$'extension',
featureNeighbourRaster, demNeighbourRaster,
args$'output-folder',
args$'dem-calc-size',
args$'exp-rast')
indic <- merge(indic, featStat, by = "trFeatId")
indicFeatCountPost <- nrow(indic)
}
# runs terrain indices
if (args$'terrain') {
ktpiFunction <- "terrain"
featTerr <- terrainIndices(folderPath, args$'tile-col', args$'tile-row', args$'extension',
featureNeighbourRaster,
demNeighbourRaster,
args$'output-folder',
args$'dem-calc-size',
args$'exp-rast')
indic <- merge(indic, featTerr, by = "trFeatId")
indicFeatCountPost <- nrow(indic)
}
# runs kernel topographic position indices
if (args$'ktpi') {
ktpiFunction <- "ktpi"
featKtpi <- ktpiIndices(folderPath, args$'tile-col', args$'tile-row', args$'extension',
featureNeighbourRaster,
demNeighbourRaster,
args$'output-folder',
args$'dem-calc-size',
args$'kernel-size',
args$'exp-rast')
indic <- merge(indic, featKtpi, by = "trFeatId")
indicFeatCountPost <- nrow(indic)
}
# runs kernel aspect indices
if (args$'kaspSlp' | args$'kaspDir' | args$'kaspSDir' |
args$'kaspCDir' | args$'kaspSlpSDir' | args$'kaspSlpCDir' |
args$'kaspSlpEle2' | args$'kaspSlpEle2SDir' | args$'kaspSlpEle2CDir' |
args$'kaspSlpLnEle' | args$'kaspSlpLnEleSlpSDir' | args$'kaspSlpLnEleSlpCDir'
) {
if (args$'kaspSlp') { ktpiFunction <- "kaspSlp" }
if (args$'kaspDir') { ktpiFunction <- "kaspDir" }
if (args$'kaspSDir') { ktpiFunction <- "kaspSDir" }
if (args$'kaspCDir') { ktpiFunction <- "kaspCDir" }
if (args$'kaspSlpSDir') { ktpiFunction <- "kaspSlpSDir" }
if (args$'kaspSlpCDir') { ktpiFunction <- "kaspSlpCDir" }
if (args$'kaspSlpEle2') { ktpiFunction <- "kaspSlpEle2" }
if (args$'kaspSlpEle2SDir') { ktpiFunction <- "kaspSlpEle2SDir" }
if (args$'kaspSlpEle2CDir') { ktpiFunction <- "kaspSlpEle2CDir" }
if (args$'kaspSlpLnEle') { ktpiFunction <- "kaspSlpLnEle" }
if (args$'kaspSlpLnEleSlpSDir') { ktpiFunction <- "kaspSlpLnEleSlpSDir" }
if (args$'kaspSlpLnEleSlpCDir') { ktpiFunction <- "kaspSlpLnEleSlpCDir" }
featKaspi <- kaspIndices(ktpiFunction,
folderPath, args$'tile-col', args$'tile-row', args$'extension',
featureNeighbourRaster,
demNeighbourRaster,
args$'output-folder',
args$'dem-calc-size',
args$'kernel-size',
args$'exp-rast',
args$'orientation')
indic <- merge(indic, featKaspi, by = "trFeatId")
indicFeatCountPost <- nrow(indic)
}
if (indicFeatCountInit != indicFeatCountPost) stop("processed stopped: input feature count is not equivalent to output feature count")
# turn scientific notation off
options(scipen=999)
# round off numeric indic values to 6 decimal places
is.num <- sapply(indic, is.numeric)
indic[is.num] <- lapply(indic[is.num], round, 6)
# replace indic "NA" values with ""
indic[is.na(indic)] <- ""
# generates output indices table filename and path
outputTableName <- paste(args$'tile-col', args$'tile-row', ktpiFunction, args$'orientation', args$'dem-calc-size', args$'kernel-size', "indices.csv", sep = "_")
outputTablePath <- paste(args$'output-folder', outputTableName, sep = "/")
# writes out indices table file
write.table(indic, file = outputTablePath, append = FALSE, sep = ",", col.names = TRUE, row.names = FALSE, quote = FALSE)
}
if (args$'neighbours') {
source_local("lib/ktpi_util.R")
neighbours <- getNeighbours(args$'folder', args$'tile-col', args$'tile-row', args$'extension',
args$'tile-col-min', args$'tile-col-max', args$'tile-row-min', args$'tile-row-max',
args$'raster-cells', args$'raster-cell-size', args$'kernel-size', tiles)
writeLines(neighbours)
}
if (args$'clicommands') {
source_local("lib/ktpi_util.R")
createKtpiCLICommands(args$'ktpi-function', args$'folder', args$'dem-folder', args$'extension', args$'output-folder',
args$'tile-col-min', args$'tile-col-max', args$'tile-row-min', args$'tile-row-max',
args$'raster-cells', args$'raster-cell-size', args$'dem-calc-size', args$'kernel-from', args$'kernel-to',
args$'kernel-step', args$'orientation', args$'exp-rast', tilesfile, tiles)
}
if (args$'jsonmessages') {
source_local("lib/ktpi_util.R")
createKtpiJSONMessages(args$'ktpi-function', args$'folder', args$'dem-folder', args$'extension', args$'output-folder',
args$'tile-col-min', args$'tile-col-max', args$'tile-row-min', args$'tile-row-max',
args$'raster-cells', args$'raster-cell-size', args$'dem-calc-size', args$'kernel-from', args$'kernel-to',
args$'kernel-step', args$'orientation', args$'exp-rast', tilesfile, tiles)
}