/
tb_bluedjinn.R
1069 lines (930 loc) · 39.4 KB
/
tb_bluedjinn.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
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#' Read a GPX file.
#' @description Read a GPX file. By default, it reads all possible GPX layers, and only returns shapes for layers that have any features.
#' if the layer has any features a sp object is returned.
#'
#' @param file a GPX filename (including directory)
#' @param layers vector of GPX layers. Possible options are \code{"waypoints"}, \code{"tracks"}, \code{"routes"}, \code{"track_points"}, \code{"route_points"}. By dedault, all those layers are read.
#' @keywords internal
#' @return a sp* object containing the gpx track objects
#' @export read_gpx
#' @note adapted from \code{\link[tmaptools]{read_GPX}}
#'
read_gpx <- function(file,
layers=c("waypoints", "tracks", "routes", "track_points", "route_points")
) {
if (!all(layers %in% c("waypoints", "tracks", "routes", "track_points", "route_points"))) stop("Incorrect layer(s)", call. = FALSE)
# check if features exist per layer
suppressWarnings(hasF <- sapply(layers, function(file,l) {
rgdal::ogrInfo(dsn = file, layer=l)$have_features
}))
if (!any(hasF)) stop("None of the layer(s) has any features.", call. = FALSE)
res <- lapply(layers[hasF], function(l) {
rgdal::readOGR(dsn = file, layer=l, verbose=FALSE)
})
names(res) <- layers[hasF]
if (sum(hasF)==1) {
res[[1]]
} else {
res
}
}
#' Read, convert and write xyz DEM/DSM Data as typically provided by the Authorities to GeoTiff.
#'
#' @description
#' Read xyz data, generate a GeoTiff file and optionally returns a raster* object.
#'
#' @param xyzFN ASCII tect file with xyz values
#' @param epsgCode "25832"
#' @param returnRaster logical. return as raster if TRUE
#' @return object of class raster*
#'
#' @export xyz2tif
#' @examples
#'\dontrun{
#' ##- libraries
#' require(uavRst)
#' setwd(tempdir())
#' ##- get typical xyz DEM data in this case from the Bavarian authority
#' utils::download.file("http://www.ldbv.bayern.de/file/zip/10430/DGM_1_ascii.zip",
#' "testdata.zip")
#' file<- unzip("testdata.zip",list = TRUE)$Name[2]
#' unzip("testdata.zip",files = file, overwrite = TRUE)
#' ##- show structure
#' head(read.table(file))
#' ##- create tiff file
#' ##- NOTE for creating a geotiff you have to provide the correct EPSG code from the meta data
#' xyz2tif(file,epsgCode = "31468")
#'
#' ##- visualize it
#' raster::plot(raster::raster(file))
#' }
#'
xyz2tif <- function(xyzFN=NULL,
epsgCode ="25832",
returnRaster = TRUE){
# read data
xyz<-data.table::fread(xyzFN)
message("this will probably take a while...\n")
r <- raster::rasterFromXYZ(xyz,crs=sp::CRS(paste0("+init=epsg:",epsgCode)))
# write it to geotiff
raster::writeRaster(r, paste0(dirname(xyzFN),"/",tools::file_path_sans_ext(basename(xyzFN)),".tif"),overwrite=TRUE)
message("...finished\n")
}
# adjust projection of objects according to their keywords
h_raster_adjust_projection <- function(x) {
llcrs <- "+proj=longlat +datum=WGS84 +no_defs"
is.fact <- raster::is.factor(x)[1]
non_proj_waning <-
paste("supplied", class(x)[1], "has no projection information!", "\n",
"provide a correctly georeferenced data raster object or 'GDAL File")
if (is.fact) {
x <- raster::projectRaster(
x, raster::projectExtent(x, crs = sp::CRS(llcrs)),
method = "ngb")
x <- raster::as.factor(x)
} else {
x <- raster::projectRaster(
x, raster::projectExtent(x, crs = sp::CRS(llcrs)),
method = "bilinear")
}
return(x)
}
# Check projection of objects according to their keywords
h_comp_ll_proj4 <- function(x) {
proj <- datum <- nodefs <- "FALSE"
allWGS84 <- as.vector(c("+init=epsg:4326", "+proj=longlat", "+datum=WGS84", "+no_defs", "+ellps=WGS84", "+towgs84=0,0,0"))
s <- as.vector(strsplit(x," "))
for (i in seq(1:length(s[[1]]))) {
if (s[[1]][i] == "+init=epsg:4326") {
proj <- datum <- nodefs <- "TRUE"
}
if (s[[1]][i] == "+proj=longlat") {
proj <- "TRUE"
}
if (s[[1]][i] == "+no_defs") {
nodefs <- "TRUE"
}
if (s[[1]][i] == "+datum=WGS84") {
datum <- "TRUE"
}
}
if (proj == "TRUE" & nodefs == "TRUE" & datum == "TRUE") {
ret <- TRUE
} else {
ret = FALSE
}
return(ret)
}
#' applies a line to a raster and returns the position of the maximum value.
#' @description applies a line to a raster and returns the position of the maximum value
#' @param dem raster object
#' @param line sp object
#' @return returns the coordinate of the maximum value of any underlying raster object along a line of interest
#' @keywords internal
#' @export
#'
line_maxpos <- function(dem,line){
mask <- dem
raster::values(mask) <- NA
#...update it with the altitude information of the flightline
mask <- raster::rasterize(line,mask)
mask2 <- mask*dem
# and find the position of the max altitude
idx = raster::which.max(mask2)
maxPos = raster::xyFromCell(mask2,idx)
return(maxPos)
}
#' extract for all polygons the position of the maximum value of the applied raster(s).
#' @description
#' extract for all polygons the position of the maximum value
#' @param fileName path and name of a GDAL raster file
#' @param layerName layer name of shape file
#' @param polySplit split polygon in single file, default is TRUE
#' extract for all polygons the position of the maximum value
#' @param cores number of cores used
#' @return returns a list of coordinates of the position of the maximum value and geometric centre of any corresponding raster object inside a given polygon
#' @keywords internal
#' @export poly_maxpos
#'
poly_maxpos <- function(fileName,layerName, polySplit=TRUE, cores=1){
if (!exists("path_run")) path_tmp = tempdir()
# read raster input data
if (polySplit) {system(paste0("rm -rf ",paste0(path_tmp,"split")))}
dem <- raster::raster(fileName)
fn <- spatial.tools::create_blank_raster(reference_raster=dem,filename = paste0(path_tmp,basename(layerName),"raw"))
mask <- raster::raster(fn)
maskx <- velox::velox(mask)
# chmx <- velox::velox(dem)
# read vector input data the sf way
sf_dcs <- sf::st_read(paste0(layerName,".shp"),quiet = TRUE)
dcs <- methods::as(sf_dcs, "Spatial")
# retrieve unique NAME
ids <- unique(dcs@data$NAME)
if (polySplit) {
message(" split polygons...\n")
message(" analyze",length(ids) ,"polygons\n")
message(" calculaton time is approx.: ",floor(length(ids)/180)," min\n")
dir.create(paste0(path_tmp,"split"),recursive=TRUE)
# split polygon with respect to the NAME attribute
parallel::mclapply(ids,function(x){
rn <- as.character(x)
gdalUtils::ogr2ogr(src_datasource_name = paste0(layerName,".shp"),
dst_datasource_name = paste0(path_tmp,"split/",basename(layerName),"_",rn,".shp"),
where = paste0("NAME='",rn,"'")
, nln = rn)
},
mc.cores = cores)
}
# parallel retrival of maxpos
message(" max height coords search...\n")
ret_max_pos <- parallel::mclapply(ids,function(x) {
rn <- as.character(x)
# create temp folder and assign it to raster
dir.create(paste0(path_tmp,rn),recursive=TRUE)
raster::rasterOptions(tmpdir=paste0(path_tmp,rn))
# read single polygon sf is even in this construct times faster
sf_shp <- sf::st_read(paste0(path_tmp,"split/",basename(layerName),"_",rn,".shp"),quiet = TRUE)
shp <- methods::as(sf_shp, "Spatial")
# reclass VALUE to 1
shp@data$VALUE <-1
# rasterize mask
maskx$rasterize(shp,field = "VALUE",band = 1)
# re-convert to raster format
m1 <- maskx$as.RasterLayer(band=1)
# get maxpos of crown area
m1 <-raster::crop(m1,c(sp::bbox(shp)[1],sp::bbox(shp)[3],sp::bbox(shp)[2],sp::bbox(shp)[4]))
d1 <-raster::crop(m1,c(sp::bbox(shp)[1],sp::bbox(shp)[3],sp::bbox(shp)[2],sp::bbox(shp)[4]))
max_pos <- raster::xyFromCell(d1,which.max(m1 * dem))
# write it to a df
df <- data.frame(x = max_pos[1], y = max_pos[2], id = rn)
# get rid of temp raster files
system(paste0("rm -rf ",paste0(path_tmp,rn)))
return(df)}, mc.cores = parallel::detectCores()-1
)
# create a spatial point data frame
max_pos <- as.data.frame(do.call("rbind", ret_max_pos))
sp::coordinates(max_pos) <- ~x+y
sp::proj4string(max_pos) <- as.character(dem@crs)
max_pos@data$id<- as.numeric(max_pos@data$id)
# create seeds file for rasterizing max_pos
# re-convert to raster format
fn <- spatial.tools::create_blank_raster(reference_raster=dem,filename = paste0(path_tmp,basename(layerName),"raw"))
mask <- raster::raster(fn)
seeds <- raster::rasterize(max_pos,mask,field="id")
seeds[seeds >= 0] <- 1
#raster::writeRaster(seeds,file.path(R.utils::getAbsolutePath(path_run),"seeds.tif"),overwrite=TRUE)
return(list(seeds,max_pos))
}
#' converts GRASS raster to Geotiff
#' @description converts GRASS raster to Geotiff
#' @param runDir path of working directory
#' @param layerName name GRASS raster
#' @param returnRaster return GRASS raster as an R raster object, default = FALSE
#' @return if returnRaster the converted GRASS raster is returned as an raster* object
#' @keywords internal
#'
grass2tif <- function(runDir = NULL, layerName = NULL, returnRaster = FALSE) {
link2GI::linkGRASS7()
rgrass7::execGRASS("r.out.gdal",
flags = c("c","overwrite","quiet"),
createopt = "TFW=YES,COMPRESS=LZW",
input = layerName,
output = paste0(runDir,"/",layerName,".tif")
)
if (returnRaster) return(raster::raster(paste0(runDir,"/",layerName,".tif")))
}
#' converts OGR to GRASS vector
#' @description converts OGR to GRASS vector
#' @param runDir path of working directory
#' @param layeName name GRASS raster
#' @keywords internal
shape2grass <- function(runDir = NULL, layerName = NULL) {
# import point locations to GRASS
rgrass7::execGRASS('v.in.ogr',
flags = c('o',"overwrite","quiet"),
input = paste0(layerName,".shp"),
output = layerName
)
}
#' converts GRASS vector to shape file
#' @description converts GRASS vector to shape file
#' @param runDir path of working directory
#' @param layerName name GRASS raster
#' @keywords internal
grass2shape <- function(runDir = NULL, layerName = NULL){
rgrass7::execGRASS("v.out.ogr",
flags = c("overwrite","quiet"),
input = layerName,
type = "line",
output = paste0(layerName,".shp")
)
}
# multiply two raster ----
funMultiply <- function(x)
{
# Note that x is received by the function as a 3-d array:
band1 <- x[,,1]
band2 <- x[,,2]
result <- band1*band2
# The output of the function should also be a 3-d array,
# even if it is a single band:
result <- array(result,dim=c(dim(x)[1],dim(x)[2],1))
return(result)
}
# returns maxpos of a raster ----
funWhichmax <- function(mask,value) {
raster::xyFromCell(value,which.max(mask * value))
}
#' Calculate descriptive statistics of raster as segmented by polygons
#'
#'@description
#' calculate statitiscs of polygon based raster extraction. Returns a spatialpolygon dataframe containing decriptive statistics
#'
#'@author Chris Reudenbach
#'
#'@param rasternames vector of raster* objects default is NULL
#'@param spdf spatial polygon dataframe default is NULL
#'@param count 0 1 switch
#'@param min 0 1 switch
#'@param max 0 1 switch
#'@param sum 0 1 switch
#'@param range 0 1 switch
#'@param mean 0 1 switch
#'@param var 0 1 switch
#'@param stddev 0 1 switch
#'@param quantile number of quantile
#'@param proj projection string
#'@param path_run run time folder for all kind of calculations, by default tempdir()
#'@param giLinks list of GI tools cli pathes, default is NULL
#'@param parallel run it parallel default is 1
#'@return data frame containing the descriptive statistic of all corresponding raster* objects for each given polygon
#'
#'@export poly_stat
#'@examples
#'\dontrun{
#' # required packages
#' require(uavRst)
#' require(link2GI)
#'
#' # create and check the links to the GI software
#' giLinks<-uavRst::linkGI(linkItems = c("saga","gdal"))
#' if (giLinks$saga$exist) {
#'
#' # get the rgb image, chm and training data
#' url <- "https://github.com/gisma/gismaData/raw/master/uavRst/data/tutorial_data.zip"
#' utils::download.file(url, file.path(tempdir(),"tutorial_data.zip"))
#' unzip(zipfile = file.path(tempdir(),"tutorial_data.zip"), exdir = tempdir())
#'
#' polyStat <- poly_stat("chm_3-3_train1",
#' spdf = "rgb_3-3_train1.shp",
#' giLinks=giLinks)
#'
#' raster::plot(polyStat)
#'}
#' ##+}
#'
poly_stat <- function(rasternames = NULL,
spdf = NULL,
count = 1,
min = 1,
max = 1,
sum = 1,
range = 1,
mean = 1,
var = 1,
stddev = 1,
quantile = 10,
parallel = 1,
proj = "+proj=utm +zone=30 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs",
path_run = tempdir(),
giLinks =NULL ) {
#message(":: run statistics...\n")
# calculate chm statistics for each crown
if (is.null(giLinks)){
giLinks <- list()
giLinks$saga <- link2GI::linkSAGA()
giLinks$gdal <- link2GI::linkGDAL()
}
gdal <- giLinks$gdal
saga <- giLinks$saga
sagaCmd<-giLinks$saga$sagaCmd
if (file.exists(file.path(R.utils::getAbsolutePath(path_run),paste0(rasternames[1],".tif"))))
proj<- raster::crs(raster::raster(file.path(R.utils::getAbsolutePath(path_run),paste0(rasternames[1],".tif"))))
# else if (file.exists(paste0(path_run,rasternames[1],".sdat")))
# proj<- sp::CRS(raster::raster(paste0(path_run,rasternames[1],".sdat")))
if (class(rasternames[1]) %in% c("RasterLayer","RasterStack", "RasterBrick"))
for (item in rasternames) {raster2sdat(rastobj = item ,path_run = R.utils::getAbsolutePath(path_run))}
else if (class(rasternames[1]) %in% c("character") & file.exists(file.path(R.utils::getAbsolutePath(path_run),paste0(rasternames[1],".tif"))))
for (item in rasternames) {
system(paste0(gdal$path,'gdalwarp -overwrite -q ',
"-of 'SAGA' ",
file.path(R.utils::getAbsolutePath(path_run),"dtm0.tif"),' ',
dstfile = file.path(R.utils::getAbsolutePath(path_run),"dtm.tif")
)
)
# gdalUtils::gdalwarp(file.path(R.utils::getAbsolutePath(path_run),paste0(item,".tif")),
# file.path(R.utils::getAbsolutePath(path_run),paste0(item,".sdat")),
# overwrite = TRUE,
# of = 'SAGA',
# verbose = FALSE)
}
if (class(spdf) %in% c("SpatialPolygonsDataFrame") ) {
rgdal::writeOGR(obj = spdf,
layer = "spdf",
driver = "ESRI Shapefile",
dsn = file.path(R.utils::getAbsolutePath(path_run)),
overwrite_layer = TRUE)
spdfshp<-"spdf.shp"
} else spdfshp <-spdf
for (i in seq(1:length(rasternames))) {
message(":: calculate ",rasternames[i], " statistics\n")
#saga_cmd shapes_grid 2 -GRIDS=/tmp/RtmpK0j1RP/run/rgb_3-3_train1.sgrd -POLYGONS=/tmp/RtmpK0j1RP/run/rgb_3-3_train1.shp -NAMING=1 -METHOD=2 -PARALLELIZED=1 -RESULT=/tmp/RtmpK0j1RP/run/rgb_3-3_train1.shp -COUNT=1 -MIN=1 -MAX=1 -RANGE=1 -SUM=1 -MEAN=1 -VAR=1 -STDDEV=1 -QUANTILE=0
# saga <- giLinks$saga
# sagaCmd<-saga$sagaCmd
# invisible(env<-RSAGA::rsaga.env(path = saga$sagaPath))
# RSAGA::rsaga.geoprocessor(lib = "shapes_grid", module = 2,
# param = list(GRIDS = file.path(R.utils::getAbsolutePath(path_run),paste0(x[i],".sgrd")),
# POLYGONS = file.path(R.utils::getAbsolutePath(path_run),spdf),
# NAMING = 1,
# METHOD = 2,
# COUNT = count,
# MIN = min,
# MAX = max,
# SUM = sum,
# RANGE = range,
# MEAN = mean,
# VAR = var,
# STDDEV = stddev,
# QUANTILE = quantile,
# PARALLELIZED = parallel,
# RESULT = file.path(R.utils::getAbsolutePath(path_run),basename(x[i]),"Stat.shp")
# ),
# show.output.on.console = FALSE,invisible = TRUE,
# env = env)
# tmp<-shp2so(file.path(R.utils::getAbsolutePath(paste0(path_run,"/",basename(spdf)))))
# tmp <- rgdal::readOGR(dsn = file.path(R.utils::getAbsolutePath(path_run)),
# layer = tools::file_path_sans_ext(basename(spdf)),
# verbose = FALSE)
# rgdal::writeOGR(obj = tmp,
# layer = tools::file_path_sans_ext(basename(spdf)),
# driver = "ESRI Shapefile",
# dsn = file.path(R.utils::getAbsolutePath(path_run)),
# overwrite_layer = TRUE,verbose = FALSE)
#
ret <- system(paste0(sagaCmd, " shapes_grid 2 ",
" -GRIDS ",file.path(R.utils::getAbsolutePath(path_run),paste0(rasternames[i],".sgrd")),
" -POLYGONS ",file.path(R.utils::getAbsolutePath(path_run),spdfshp),
" -NAMING 0",
" -METHOD 2",
" -COUNT ", count,
" -MIN ", min,
" -MAX ", max,
" -SUM ",sum,
" -RANGE ",range,
" -MEAN ", mean,
" -VAR ",var,
" -STDDEV ",stddev,
" -QUANTILE ",quantile,
" -PARALLELIZED ",parallel,
" -RESULT ",file.path(R.utils::getAbsolutePath(path_run),paste0(basename(rasternames[i]),"Stat.shp"))),
intern = TRUE)
stat1<-shp2so(file.path(R.utils::getAbsolutePath(paste0(path_run,"/",paste0(basename(rasternames[i]),"Stat")))))
# stat1 <- rgdal::readOGR(dsn = file.path(R.utils::getAbsolutePath(path_run)),
# layer = paste0(basename(rasternames[i]),"Stat"),
# verbose = TRUE)
names(stat1) <- gsub(names(stat1),pattern = "\\.",replacement = "")
#TODO
if (nchar(rasternames[i]) > 6 ) rasternames[i] <-substr(rasternames[i],1,6)
if (i<10) tmp_names <- gsub(names(stat1),pattern = paste0("G0",i),replacement = rasternames[i])
else if (i<10) tmp_names <- gsub(names(stat1),pattern = paste0("G",i),replacement = rasternames[i])
names(stat1)<-substr(tmp_names , 1,10)
#raster::shapefile(stat1,file.path(R.utils::getAbsolutePath(path_run),"polystat"),overwrite = TRUE)
if (i == 1) {
stat <- stat1
} else {
stat@data <- cbind(stat@data,stat1@data[4:length(names(stat1))])
#stat <- stat1
}
}
sp::proj4string(stat)<-proj
so2shp(stat,paste0(basename(rasternames[i]),"Stat"),file.path(R.utils::getAbsolutePath(paste0(path_run))))
return(stat)
}
# fill holes
fillGaps<- function (folder,layer){
message(":: fill data gaps using gdal_fillnodata... \n")
# fill data holes
if (Sys.info()["sysname"] == "Windows"){
ret <- system2(command = "gdal_fillnodata.py ",args =
paste0(folder,"/",layer,".tif ",
folder,"/",layer,".tif "))
} else {
ret <- system(paste0("gdal_fillnodata.py ", folder,layer,".tif ",
folder,layer,".tif "),intern = TRUE)
}
# write filled data back to GRASS
rgrass7::execGRASS('r.in.gdal', flags=c('o',"overwrite"), input=paste0(folder,"/",layer,".tif"), output=layer, band=1)
}
# creates names and ranges from a simple list for zrange cuts ---
makenames<-function(zr ) {
class<-list()
zrange<-list()
for ( i in 1:(length(zr[[1]]))) {
if (i == length(zr[[1]])){
class[[i]]<-c(paste0('class',zr[[1]][1],zr[[1]][length(zr[[1]])]))
zrange[[i]]<-c(zr[[1]][1], zr[[1]][length(zr[[1]])])
} else {
class[[i]]<-c(paste0('class',zr[[1]][i],zr[[1]][i+1]))
zrange[[i]]<-c(zr[[1]][i],zr[[1]][i+1])
}
}
return(list(unlist(class),zrange))
}
# extract pixel values according to an overlay and response name ---
extractTrainPixelValues <- function(imgStack=NULL,trainData=NULL,responseCol=NULL){
#extract training Area pixel values
dfTpv = data.frame(matrix(vector(), nrow = 0, ncol = length(names(imgStack)) + 1))
for (i in 1:length(unique(trainData[[responseCol]]))){
category <- unique(trainData[[responseCol]])[i]
message("\n extracting cat: ",levels(category)[i]," no: ",i," of: ",length(unique(trainData[[responseCol]])))
categorymap <- trainData[trainData[[responseCol]] == category,]
dataSet <- raster::extract(imgStack, categorymap)
dataSet <- lapply(dataSet, function(x){cbind(x, class = as.numeric(rep(category, nrow(x))))})
df <- do.call("rbind", dataSet)
dfTpv <- rbind(dfTpv, df)
}
names(dfTpv)<-gsub(names(dfTpv),pattern = "\\.",replacement = "_")
return(dfTpv)
}
# calculate the mode ---
Mode <- function(x) {
ux <- unique(x)
ux[which.max(tabulate(match(x, ux)))]
}
### getPopupStyle creates popup style -------
getPopupStyle <- function() {
# htmlTemplate <- paste(
# "<html>",
# "<head>",
# "<style>",
# "#popup",
# "{font-family: Arial, Helvetica, sans-serif;width: 20%;border-collapse: collapse;}",
# "#popup td {font-size: 1em;border: 0px solid #85ADFF;padding: 3px 20px 3px 3px;}",
# "#popup tr.alt td {color: #000000;background-color: #F0F5FF;}",
# "#popup tr.coord td {color: #000000;background-color: #A8E6A8;}",
# "div.scrollableContainer {max-height: 200px;max-width: 100%;overflow-y: auto;overflow-x: auto;margin: 0px;background: #D1E0FF;}",
# "</style>",
# "</head>",
# "<body>",
# "<div class='scrollableContainer'>",
# "<table class='popup scrollable'>",
# "<table id='popup'>")
# return(htmlTemplate)
fl <- function() system.file("templates/popup.brew", package = "mapview")
pop <- readLines(fl)
end <- grep("<%=pop%>", pop)
return(paste(pop[1:(end-2)], collapse = ""))
}
#' Split multiband image to single band SAGA files
#' @description Split multiband image to single band SAGA files. If a reference file is given, it performs a resample if necessary to avoid the numerical noise problem of SAGA extent.
#' @param fn character. filename
#' @param bandname character. list of bandnames c("red","green","blue")
#' @param startBand numerical. first band to export
#' @param startBand numerial. last band to export
#' @param refFn character. reference image for resampling
#' @param returnRaster logical. return as raster
#' @param gdalLinks list of gdal installations
#' @name split2SAGA
#' @return if returnRaster the splitted SAGA raster is returned as an raster* object
#' @keywords internal
#'@export
split2SAGA<-function(fn=NULL,
bandName=NULL,
startBand= 1,
endBand =3,
refFn=NULL,
gdalLinks=NULL,
returnRaster=FALSE){
if (is.null(gdalLinks)) gdal<- link2GI::linkGDAL()
else gdal<- gdalLinks
flist<-list()
if (!exists("path_run")) path_run = tempdir()
for (i in seq(startBand:endBand)){
outFn<-file.path(R.utils::getAbsolutePath(path_run),paste0(bandName[i],".sdat"))
#raster::writeRaster(raster::raster(fn),outFn,overwrite = TRUE,NAflag = 0,process="text")
if (!is.null(refFn))
r<-raster::raster(refFn)
else
r<-raster::raster(fn)
system(paste0(gdal$path,'gdalwarp -overwrite -q ',
"-of 'SAGA' ",
"-r 'bilinear'",' ',
"-b ",as.character(i),' ',
"-a_nodata = 0",' ',
"-tr ",paste0(raster::xres(r)," ",
raster::xres(r)),' ',
"-a_srs ", as.character(r@crs),' ',
fn[[i]]@file@name,' ',
outFn
)
)
# res<-gdalUtils::gdal_translate(src_dataset = fn[[i]]@file@name,
# dst_dataset = outFn,
# tr= paste0(raster::xres(r)," ",
# raster::xres(r)),
# b = as.character(i),
# of = "SAGA",
# a_nodata = 0,
# a_srs = as.character(r@crs) )
#
r<-raster::writeRaster(raster::resample(raster::raster(outFn),raster::raster(refFn)),
filename = outFn,
NAflag = 0,
format="SAGA",
overwrite=TRUE,progress="text")
flist<-append(flist, r)
}
if (returnRaster) return(raster::stack(flist))
}
#' colorize the cat outputs
#'@description colorize the cat outputs
#'@export
#'@keywords internal
getCrayon<-function(){
head <- crayon::black $ bgGreen
err <- crayon::red $ bold
note <- crayon::blue $ bold
ok <- crayon::green $ bold
return(list(note,err,ok,head))
}
#' create name vector corresponding to the training image stack
#' create vector containing the names of the image stack claculated using \code{\link{calc_ext}}
#' @param rgbi character. codes of the RGB indices
#' @param bandNames character. band names
#' @param stat character. stat codes
#' @param morpho character. morpho codes
#' @param edge character. edge codes
#' @param rgbTrans character. rgbTrans codes
#' @param dem charater. dem codes
#' @param l_raster number of raster layer that exist
#' @keywords internal
#'
#' @export make_bandnames
make_bandnames <- function(rgbi = NA,
bandNames = NA,
stat = FALSE,
morpho = NA,
edge = NA ,
rgbTrans = NA,
dem = NA,
l_raster = NA,
pca=NA){
if (!is.na(rgbi[1])) bandNames <- append(c("red","green","blue"),rgbi)
if (!is.na(bandNames[1])) {
if(bandNames[1] == "simple"){
bandNames <- c("Energy", "Entropy", "Correlation",
"Inverse_Difference_Moment", "Inertia",
"Cluster_Shade", "Cluster_Prominence",
"Haralick_Correlation")
} else if(bandNames[1] == "advanced"){
bandNames <- c("Hara_Mean", "Hara_Variance", "Dissimilarity",
"Sum_Average",
"Sum_Variance", "Sum_Entropy",
"Difference_of_Variances",
"Difference_of_Entropies",
"IC1", "IC2")
} else if(bandNames[1] == "higher"){
bandNames <- c("Short_Run_Emphasis",
"Long_Run_Emphasis",
"Grey-Level_Nonuniformity",
"Run_Length_Nonuniformity",
"Run_Percentage",
"Low_Grey-Level_Run_Emphasis",
"High_Grey-Level_Run_Emphasis",
"Short_Run_Low_Grey-Level_Emphasis",
"Short_Run_High_Grey-Level_Emphasis",
"Long_Run_Low_Grey-Level_Emphasis",
"Long_Run_High_Grey-Level_Emphasis")
} else if(bandNames[1] == "all"){
bandNames <- c("Energy", "Entropy", "Correlation",
"Inverse_Difference_Moment", "Inertia",
"Cluster_Shade", "Cluster_Prominence",
"Haralick_Correlation",
"Hara_Mean", "Hara_Variance", "Dissimilarity",
"Sum_Average",
"Sum_Variance", "Sum_Entropy",
"Difference_of_Variances",
"Difference_of_Entropies",
"IC1", "IC2",
"Short_Run_Emphasis",
"Long_Run_Emphasis",
"Grey-Level_Nonuniformity",
"Run_Length_Nonuniformity",
"Run_Percentage",
"Low_Grey-Level_Run_Emphasis",
"High_Grey-Level_Run_Emphasis",
"Short_Run_Low_Grey-Level_Emphasis",
"Short_Run_High_Grey-Level_Emphasis",
"Long_Run_Low_Grey-Level_Emphasis",
"Long_Run_High_Grey-Level_Emphasis")
}
if (!is.na(l_raster)) bandNames<-bandNames[1:l_raster]
}
if (stat == TRUE) {
bandNames = c("Stat_Mean","Stat_Variance", "Stat_Skewness", "Stat_Kurtosis")
}
if (!is.na(dem)) {
bandNames = dem
}
if (!is.na(morpho)) {
bandNames = morpho
}
if (!is.na(edge)) {
bandNames = edge
}
if (!is.na(pca)[1]) {
bandNames = "PCA"
}
if (!is.na(rgbTrans)) {
if (rgbTrans %in% c("Gray"))
bandNames = bandNames <- c(paste0(rgbTrans,"_b1"))
# c("cielab","CMY","Gray","HCL","HSB","HSI","Log","XYZ","YUV")
else if (rgbTrans %in% c("cielab"))
bandNames = bandNames <- c(paste0(rgbTrans,"_b1"),paste0(rgbTrans,"_b2"),paste0(rgbTrans,"_b3"),paste0(rgbTrans,"_b4"))
else
bandNames = bandNames <- c(paste0(rgbTrans,"_b1"),paste0(rgbTrans,"_b2"),paste0(rgbTrans,"_b3"))
}
return(bandNames)
}
# returns the saga items from a list --
issagaitem <- function(x)
{
if (x %in% c("SLOPE","ASPECT","C_GENE","C_PROF","C_PLAN","C_TANG","C_LONG","C_CROS","C_MINI","C_MAXI","C_TOTA","C_ROTO","MTPI") ) return(TRUE) else return(FALSE)
}
# returns the gdal items from a list ---
isgdaldemitem <- function(x)
{
if (x %in% c("hillshade","slope", "aspect","TRI","TPI","Roughness")) return(TRUE) else return(FALSE)
}
# if (substr(Sys.getenv("COMPUTERNAME"),1,5) == "PCRZP") {
# gdalUtils::gdal_setInstallation(search_path = shQuote("C:/Program Files/QGIS 2.14/bin/"))
# } else {
# ## (gdalUtils) check for a valid GDAL binary installation on your system
# if (!quiet) gdalUtils::gdal_setInstallation(verbose = TRUE)
# else gdalUtils::gdal_setInstallation()
# }
#'@title Checks if running on a specified computer domain
#'@name setHomePath
#'@description Checks if the computer name belongs to a specified group i.e. aq network domain Marburg Universitys computer pools
#'@param homeDir full path to the real folder location the project
#'@param prefixPC contains an arbitrary part of the computer name. It always starts with the first letter.
#'@author CR
#'@keywords internal
#' @return string containing a directory path
#'@examples
#' \dontrun{
#' # add path
#' setHomePath("saga",prefixPC="PCRZP")
#' }
#'@export setHomePath
setHomePath<- function(homeDir="F:/MPG", prefixPC="PCRZP") {
if (!exists("GiEnv")) GiEnv <- new.env(parent=globalenv())
if (substr(Sys.getenv("COMPUTERNAME"),1,nchar(prefixPC)) == substr(prefixPC,1,nchar(prefixPC))) {
projHomeDir <- shQuote(homeDir)
return(path.expand(projHomeDir))
} else {
return(path.expand("~/edu"))
}
}
#' returns the status of saga modules
#' @param module character name of module to be checked
#' @param file character. filename to be checked
#' @param listname character. name of log list
#'@keywords internal
#'@export
fileProcStatus <- function(module=NULL,file= NULL,listname=NULL){
if (!exists("path_run")) path_run = tempdir()
assign(listname,list())
if (file.exists(file.path(R.utils::getAbsolutePath(path_run),file))) {
eval(parse(text=paste0(listname,"$",module," <- TRUE")))
return(eval(parse(text=listname)))
}
else {
eval(parse(text=paste0(listname,"$",module," <- FALSE")))
return(eval(parse(text=listname)))
message(eval(parse(text=listname)))
}
}
#' writes shapefiles from sf or sp* objects
#' @param spobj spatial object of type sp* or sf
#' @param name character. name of object to be created
#' @param path_run character. path used for runtime operations
#'@keywords internal
#'@export
so2shp<-function(spobj,
name="spobj",
path_run = tempdir()){
if(class(spobj[[1]]) %in% c("SpatialPolygons","SpatialPolygonsDataFrame","SpatialPointsDataFrame","SpatialLinesDataFrame","SpatialLines","SpatialPoints")) {
raster::shapefile(spobj,file.path(R.utils::getAbsolutePath(path_run),paste0(name,".shp")))
} else if(class(spobj)[[1]] == "sf") {
# read vector input data the sf way
sf::st_write(obj=spobj,dsn= path_run ,layer=paste0(name,".shp"),driver="ESRI Shapefile",quiet = TRUE, delete_dsn=TRUE)
}
}
#' reads shapefiles as sf or sp objects
#' @param fn filename optionally with path
#' @param type character. type of object to be created. "sp" is default can be set to "sf"
#' @param path_run character. path used for runtime operations
#' @return a sf* or sf* object from a given shape file
#'@keywords internal
#'@export
shp2so<-function(fn,
type="sp",
path_run = tempdir()){
fname<- tools::file_path_sans_ext(basename(fn))
if (nchar(dirname(fn))>1) path_run <- dirname(fn)
if(type=="sp") {
spobj<-raster::shapefile(file.path(R.utils::getAbsolutePath(path_run),paste0(fname,".shp")))
} else if(class(spobj)[[1]] == "sf") {
# read vector input data the sf way
spobj<-sf::st_read(dsn= file.path(path_run ,layer=paste0(fname,".shp")),quiet = TRUE)
}
return(spobj)
}
#' writes raster* objects to the SAGA format
#' @param rastobj raster* object
#' @param path_run character. path used for runtime operations
#'@keywords internal
raster2sdat<-function(rastobj,
path_run = tempdir()
){
if (class(rastobj) %in% c("RasterLayer")) {
raster::writeRaster(rastobj,file.path(R.utils::getAbsolutePath(path_run),paste0(names(rastobj),".sdat")),overwrite = TRUE,NAflag = 0)
} else if (class(rastobj) %in% c("RasterStack", "RasterBrick")) {
i=1
for (item in names(rastobj)){
raster::writeRaster(rastobj[[i]],
file.path(R.utils::getAbsolutePath(path_run),paste0(item,".sdat")),
overwrite = TRUE,
NAflag = 0)
i=i+1
}
}
}
#' converts SAGA raster to R raster object
#' @description converts SAGA raster to R raster object
#' @param fn filname without extension
#' @param ext extent of the raster in R notation
#' @param path_run character. path used for runtime operations
#' @param gdalLinks list of gdal installations
#' @return returned saga raster file as an raster* object
#' @keywords internal
#' @export
saga2r<- function(fn,
ext,
path_run=tempdir(),
gdalLinks=NULL ) {
if (is.null(gdalLinks)) gdal<- link2GI::linkGDAL()
else gdal<- gdalLinks
system(paste0(gdal$path,'gdalwarp -overwrite -q ',
"-of 'SAGA' ",
file.path(R.utils::getAbsolutePath(path_run),fn,".sdat"),' ',
dstfile = file.path(R.utils::getAbsolutePath(path_run),fn,".tif")
)
)
# gdalUtils::gdalwarp(paste0(path_run,fn,".sdat"),
# paste0(path_run,fn,".tif"),
# overwrite = TRUE,
# verbose = FALSE)
x<-raster::raster(paste0(path_run,fn,".tif"))
x@extent <- ext
# convert to SAGA
return(x)
}
#' converts R raster* objects to SAGA raster
#' @param x raster object
#' @param fn filname name without extension
#' @param path_run character. path used for runtime operations
#' @keywords internal
#' @export
r2saga <- function(x,fn,path_run=tempdir()) {
fname<- tools::file_path_sans_ext(basename(fn))
dirname<-dirname(fn)
if (nchar(dirname)>1) path_run<-dirname
raster::writeRaster(x,file.path(R.utils::getAbsolutePath(path_run),paste0(fn,".sdat")),bylayer=TRUE,overwrite = TRUE,NAflag = 0)
}
#' convenient function to establish all link2GI links
#' @description brute force search, find and linkl of all link2GI link functions. This is helpfull if yor system is wellsetup and the standard linkage procedure will provide the correct linkages.
#'
#' @note You may also use the full list of arguments that is made available from the \code{link2GI} package, but it is strongly recommended in this case to use directly the single linkage functions from \code{link2GI}.
#' @param links character. links
#' @param linkItems character. list of c("saga","grass7","otb","gdal")
#' @param simple logical. true make all
#' @param sagaArgs character. full string of sagaArgs
#' @param grassArgs character. grassArgs full string of grassArgs
#' @param otbArgs character. full string of otbArgs
#' @param gdalArgs character. full string of gdalArgs
#' @param quiet supress all messages default is FALSE
#' @return list of all link2GI bindings
#'@examples