1 |
7526fb1c
|
Jim Regetz
|
# 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")
|