1
|
# Calculates the monthly average clear day coverage for the clear day
|
2
|
# MODIS layer extracted from the MODIS LST files
|
3
|
#
|
4
|
# Reads a directory of .img files and groups them by month to calculate
|
5
|
# a monthly aveage. This script also applies the clear day scaling factor
|
6
|
# to the resulting average raster and writes one raster for each month
|
7
|
# of the 1-year 2000 to 2010 period (one file for each month of the
|
8
|
# 10-year group)
|
9
|
#
|
10
|
# Developed by John Donoghue
|
11
|
# Created: 20 October 2010
|
12
|
# Last Updated:
|
13
|
#
|
14
|
# For NCEAS Working Group Environment and Organisms
|
15
|
#
|
16
|
|
17
|
#setwd("/data/project/organisms/R")
|
18
|
#source("Calc Clear Day Monthly Avg.r", echo=TRUE, print.eval=TRUE)
|
19
|
|
20
|
wkDir <- "/data/project/organisms/MODIS_LST_Oregon/ClearDayGDAL/ClearDay_Original_IMG_Extracts/ByDate"
|
21
|
setwd(wkDir)
|
22
|
|
23
|
inFiles <- list.files(pattern="*h08v04.005.ClearDay.img$")
|
24
|
nFiles <- length(inFiles)
|
25
|
|
26
|
library(raster)
|
27
|
|
28
|
# build a list of the files for each month
|
29
|
janList <- list()
|
30
|
febList <- list()
|
31
|
marList <- list()
|
32
|
aprList <- list()
|
33
|
mayList <- list()
|
34
|
junList <- list()
|
35
|
julList <- list()
|
36
|
augList <- list()
|
37
|
sepList <- list()
|
38
|
octList <- list()
|
39
|
novList <- list()
|
40
|
decList <- list()
|
41
|
|
42
|
# for each file in the list
|
43
|
for (i in 1:nFiles) {
|
44
|
|
45
|
#print(inFiles[i])
|
46
|
|
47
|
# get month of file
|
48
|
theMonth <- substr(inFiles[i], 16,17)
|
49
|
if (substr(theMonth, 2,2) == ".") {theMonth <- substr(theMonth, 1, 1)}
|
50
|
|
51
|
# group by month
|
52
|
if (theMonth == "1") {
|
53
|
a <- length(janList)
|
54
|
janList[a+1] <- inFiles[i]
|
55
|
}
|
56
|
|
57
|
if (theMonth == "2") {
|
58
|
a <- length(febList)
|
59
|
febList[a+1] <- inFiles[i]
|
60
|
}
|
61
|
|
62
|
if (theMonth == "3") {
|
63
|
a <- length(marList)
|
64
|
marList[a+1] <- inFiles[i]
|
65
|
}
|
66
|
|
67
|
if (theMonth == "4") {
|
68
|
a <- length(aprList)
|
69
|
aprList[a+1] <- inFiles[i]
|
70
|
}
|
71
|
|
72
|
if (theMonth == "5") {
|
73
|
a <- length(mayList)
|
74
|
mayList[a+1] <- inFiles[i]
|
75
|
}
|
76
|
|
77
|
if (theMonth == "6") {
|
78
|
a <- length(junList)
|
79
|
junList[a+1] <- inFiles[i]
|
80
|
}
|
81
|
|
82
|
if (theMonth == "7") {
|
83
|
a <- length(julList)
|
84
|
julList[a+1] <- inFiles[i]
|
85
|
}
|
86
|
|
87
|
if (theMonth == "8") {
|
88
|
a <- length(augList)
|
89
|
augList[a+1] <- inFiles[i]
|
90
|
}
|
91
|
|
92
|
if (theMonth == "9") {
|
93
|
a <- length(sepList)
|
94
|
sepList[a+1] <- inFiles[i]
|
95
|
}
|
96
|
|
97
|
if (theMonth == "10") {
|
98
|
a <- length(octList)
|
99
|
octList[a+1] <- inFiles[i]
|
100
|
}
|
101
|
|
102
|
if (theMonth == "11") {
|
103
|
a <- length(novList)
|
104
|
novList[a+1] <- inFiles[i]
|
105
|
}
|
106
|
|
107
|
if (theMonth == "12") {
|
108
|
a <- length(decList)
|
109
|
decList[a+1] <- inFiles[i]
|
110
|
}
|
111
|
|
112
|
# end list files
|
113
|
}
|
114
|
|
115
|
print("")
|
116
|
print("PROCESSING h08v04 FILES")
|
117
|
print("")
|
118
|
|
119
|
# for each file in the list
|
120
|
for (f in 1:12) {
|
121
|
if (f == 1) {myList <- janList}
|
122
|
if (f == 2) {myList <- febList}
|
123
|
if (f == 3) {myList <- marList}
|
124
|
if (f == 4) {myList <- aprList}
|
125
|
if (f == 5) {myList <- mayList}
|
126
|
if (f == 6) {myList <- junList}
|
127
|
if (f == 7) {myList <- julList}
|
128
|
if (f == 8) {myList <- augList}
|
129
|
if (f == 9) {myList <- sepList}
|
130
|
if (f == 10) {myList <- octList}
|
131
|
if (f == 11) {myList <- novList}
|
132
|
if (f == 12) {myList <- decList}
|
133
|
|
134
|
print ("")
|
135
|
print ("")
|
136
|
|
137
|
filestoStack <- raster()
|
138
|
for (iCtr in 1 : length(myList))
|
139
|
{
|
140
|
mFile <- myList[iCtr]
|
141
|
print(sprintf("..Reading file: %s",mFile))
|
142
|
myFile <- paste(wkDir, "/", mFile, sep="")
|
143
|
inputRaster <- raster()
|
144
|
inputRaster <- mFile
|
145
|
filestoStack <- c(filestoStack, inputRaster)
|
146
|
}
|
147
|
|
148
|
# now calculate the mean cell values
|
149
|
print("....Stacking Grids")
|
150
|
myStack <- stack(filestoStack)
|
151
|
|
152
|
# Check if file exists
|
153
|
outFile <- paste("Month_", f ,"_ClearDay_Average",".img",sep="")
|
154
|
if (file.exists(outFile) == FALSE) {
|
155
|
print("....Averaging Monthly Values")
|
156
|
meanRaster <- mean(myStack, na.rm=TRUE)
|
157
|
print("....Applying Scale Factor")
|
158
|
scaledRaster <- meanRaster * 0.0005 # apply scaling factor to raster
|
159
|
print(sprintf("......Writing file: %s",outFile))
|
160
|
writeRaster(scaledRaster,filename=outFile,format="HFA",overwrite=TRUE)
|
161
|
}
|
162
|
}
|
163
|
|
164
|
# finished
|
165
|
print("Finished")
|
166
|
|