Newer
Older
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
#' Unify extent of a rasterstack
max_extent <- function(rl){
# given list of rasters
# returns union of extent
xmin=min(sapply(rl,function(x){extent(x)@xmin}))
xmax=max(sapply(rl,function(x){extent(x)@xmax}))
ymin=min(sapply(rl,function(x){extent(x)@ymin}))
ymax=max(sapply(rl,function(x){extent(x)@ymax}))
extent(c(xmin,xmax,ymin,ymax))
}
# Saving function
writeGeoTiff <- function(file, fname,dt = "INT2S"){
require(assertthat)
if(!has_extension(fname,"tif")) fname <- paste0(fname,".tif")
writeRaster(file,fname,
format='GTiff',
datatype = dt,
NAflag = -9999,
options=c("COMPRESS=DEFLATE","PREDICTOR=2","ZLEVEL=9"),
overwrite= TRUE )
}
emptyraster <- function(x, ...) { # add name, filename,
emptyraster <- raster(nrows=nrow(x), ncols=ncol(x),
crs=x@crs,
ext=extent(x), ...)
return(emptyraster)
}
#' Align raster data by bringing it in the same geometry and extend.
alignRasters <- function(data, template, method = "bilinear",func = mean,cl = T,expand = TRUE){
lib <- c("raster")
sapply(lib, function(...) stopifnot(require(..., character.only = T)))
if(cl) beginCluster(parallel::detectCores()-1)
if(projection(data) == projection(template)){
data <- raster::crop(data, template, snap = "near")
if(class(template) == "RasterLayer"){
if(data@ncols / template@ncols >= 2){
factor <- floor(data@ncols/template@ncols)
data <- raster::aggregate(data, fact = factor, fun = func,na.rm = TRUE,
expand=expand)
}
data <- raster::resample(data, template, method = method)
}
} else {
data <- projectRaster(data, template, method = method)
}
if(cl) endCluster()
return(data)
}
split_raster <- function(infile,outpath,parts=3) {
## parts = division applied to each side of raster, i.e. parts = 2 gives 4 tiles, 3 gives 9, etc.
lib <- c("gdalUtils","parallel","reshape2")
sapply(lib, function(...) stopifnot(require(..., character.only = T)))
# Check if infile exists
stopifnot(file.exists(infile))
filename <- strsplit(basename(infile),".tif")[[1]]
# Get the dimensions of the file
dims <- as.numeric(
strsplit(gsub('Size is|\\s+', '', grep('Size is', gdalinfo(infile), value=TRUE)),
',')[[1]]
)
# Generate window frames
xy <- list()
# t is nr. of iterations per side
t <- parts - 1
for (i in 0:t) {
for (j in 0:t) {
# [-srcwin xoff yoff xsize ysize] src_dataset dst_dataset
srcwin <- paste(i * dims[1]/parts, j * dims[2]/parts, dims[1]/parts, dims[2]/parts)
xy[[paste0(i,"_",j)]] <- data.frame(infile = infile,srcwin = srcwin, file = paste0(outpath,"/",filename,"_",i,"_",j,".tif"))
}
}
df <- melt(xy)
# Then process per src_win
cat("Start splitting: ",filename)
# Create a function to split the raster using gdalUtils::gdal_translate
split <- function(input, outfile, srcwin) {
gdal_translate(input, outfile, srcwin=srcwin)
}
# Make a copy for export
df_org <- df
# Kick out files already existing
df <- df[which(!file.exists(as.character(df$file))),]
# Make cluster
cl <- makeCluster(detectCores()-1)
clusterExport(cl, c('split', 'df'))
clusterEvalQ(cl,library(gdalUtils))
system.time({
parLapply(cl, seq_len(nrow(df)), function(i) {
split(df$infile[i], df$file[i], df$srcwin[i])
})
})
stopCluster(cl)
cat("\n")
cat("Done")
return(df_org)
}