-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathstartupscript1b.R
94 lines (77 loc) · 2.98 KB
/
startupscript1b.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
# Author: B.P. Ottow
# parameters
setwd("/media/boukepieter/schijfje_ottow/thesis/workspace1b")
channelLength = 500
minimumHead = 25
minimumDebiet = 0.1
minimumPotential = 100000
# files
ETdir <- "input/ET"
Pdir <- "input/P"
DEMfile <- "input/DEMs.tif"
source("model1b.R")
HydroPowerMonthly(DEMfile, Pdir, ETdir, minimumPotential = minimumPotential, minimumDebiet = minimumDebiet,
minimumHead = minimumHead, channelLength = channelLength)
####################################-Preprosessing-##################################################
library(raster)
install.packages("gdalUtils")
library(gdalUtils)
utm51 = CRS('+proj=utm +zone=51 +ellps=WGS84 +datum=WGS84 +units=m +no_defs')
wgs84 = CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
files <- list.files("Data/MOD16")
dir.create("monthly")
hdf4_dataset <- system.file(paste("Data/MOD16/", files[1], sep=""), package="gdalUtils")
for (i in 1:length(files)){
gdal_translate(paste("Data/MOD16/", files[i], sep=""),
paste("Data/ET/", substring(files[i],1,42), "tif",sep=""), sd_index=1)
print(i)
}
files <- list.files("Data/ET")
days <- c(31,28,31,30,31,30,31,31,30,31,30,31)
for (i in 1:12){
mod2010 <- raster(paste("Data/ET/", files[i], sep=""))
mod2011 <- raster(paste("Data/ET/", files[i+12], sep=""))
mod2012 <- raster(paste("Data/ET/", files[i+24], sep=""))
mod2013 <- raster(paste("Data/ET/", files[i+36], sep=""))
#mod2010s <- crop(mod2010, arealarge)
mod2010[mod2010>3000] <- NA
#mod2011s <- crop(mod2011, arealarge)
mod2011[mod2011>3000] <- NA
#mod2012s <- crop(mod2012, arealarge)
mod2012[mod2012>3000] <- NA
#mod2013s <- crop(mod2013, arealarge)
mod2013[mod2013>3000] <- NA
mod16 <- (mod2010 + mod2011 + mod2012 + mod2013) / 4 / days[i]
mod16utm = projectRaster(mod16, crs=utm51)
ET <- crop(mod16utm, projectarea)
writeRaster(ET, filename=paste("input/ET", i, ".tif", sep=""), format="GTiff", overwrite=TRUE)
print(i)
}
# TRMM
setwd("/media/boukepieter/schijfje_ottow/thesis/data/TRMM/3-hourly-tiff")
days <- c(31,28,31,30,31,30,31,31,30,31,30,31)
files <- list.files("monthly")
for (i in 1:12){
setwd("/media/boukepieter/schijfje_ottow/thesis/data/TRMM/3-hourly-tiff")
trmm2010 <- raster(paste("monthly/", files[i], sep=""))
trmm2011 <- raster(paste("monthly/", files[i+12], sep=""))
trmm2012 <- raster(paste("monthly/", files[i+24], sep=""))
trmm2013 <- raster(paste("monthly/", files[i+36], sep=""))
trmm <- (trmm2010 + trmm2011 + trmm2012 + trmm2013) / 40 * 3 / days[i]
projection(trmm) <- wgs84
trmmutm = projectRaster(trmm, crs=utm51)
P <- crop(trmmutm, projectarea)
setwd("/media/boukepieter/schijfje_ottow/thesis/workspace1b")
writeRaster(P, filename=paste("input/P", i, ".tif", sep=""), format="GTiff", overwrite=TRUE)
print(i)
}
DEM <- raster("Data/DEM.tif")
test <- raster("input/P1.tif")
plot(test, main="August", col=rainbow(12))
for (i in 2:12) {
m <- raster(paste("input/P", i, ".tif", sep=""))
test <- test + m
}
test <- test / 12
plot(test)
plot(P)