1
|
# Calculates the daily 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 day to calculate
|
5
|
# a daily aveage. This script also applies the clear day scaling factor
|
6
|
# to the resulting average raster and writes one raster for each day
|
7
|
# of the 1-year 2000 to 2010 period (one file for each 10-year group)
|
8
|
#
|
9
|
# Developed by John Donoghue
|
10
|
# Created: 16 October 2010
|
11
|
# Last Updated:
|
12
|
#
|
13
|
# For NCEAS Working Group Environment and Organisms
|
14
|
#
|
15
|
|
16
|
#setwd("/data/project/organisms/R")
|
17
|
#source("Calc Clear Day Daily Avg.r", echo=TRUE, print.eval=TRUE)
|
18
|
|
19
|
setwd("/data/project/organisms/MODIS_LST_Oregon/ClearDayGDAL")
|
20
|
inFiles <- list.files(pattern="*.img$")
|
21
|
|
22
|
library(raster)
|
23
|
|
24
|
print("")
|
25
|
print("Processing h08v04 tiles")
|
26
|
print("")
|
27
|
|
28
|
# process h08v04 tiles
|
29
|
# Create day averages for each day in the year
|
30
|
for(i in 1:366){
|
31
|
if (i < 9) {f <- "00"}
|
32
|
if ((i > 9) & (i <99 )) {f <- "0"}
|
33
|
if (i > 99) {f <- ""}
|
34
|
fpart <- paste(toString(f), toString(i), sep="")
|
35
|
|
36
|
# get a list of all files for that day of the year
|
37
|
print("")
|
38
|
print(sprintf("Reading day: %s", toString(fpart)))
|
39
|
myFiles <- grep(fpart,substr(inFiles, 14,16))
|
40
|
nFiles <- length(myFiles)
|
41
|
|
42
|
# for each file in the list
|
43
|
filestoStack <- list()
|
44
|
for (iCtr in 1 : nFiles)
|
45
|
{
|
46
|
mFile <- myFiles[iCtr]
|
47
|
if (substr(inFiles[mFile],18,23) == "h08v04") {
|
48
|
print(sprintf("..Reading file: %s",inFiles[mFile]))
|
49
|
inputRaster <- raster()
|
50
|
inputRaster <- raster(inFiles[mFile])
|
51
|
filestoStack <- c(filestoStack, inFiles[mFile])
|
52
|
}
|
53
|
}
|
54
|
# now calculate the mean cell values
|
55
|
myStack <- stack(filestoStack)
|
56
|
print(sprintf("....Averaging Daily Values"))
|
57
|
|
58
|
# Check if file exists
|
59
|
outFile <- paste("Day_",toString(i),"_Average_Scaled_h08v04",".img",sep="")
|
60
|
if (file.exists(outFile) == FALSE) {
|
61
|
meanRaster <- mean(myStack, na.rm=TRUE)
|
62
|
scaledRaster <- meanRaster * 0.0005 # apply scaling factor to raster
|
63
|
print(sprintf("......Writing file: %s",outFile))
|
64
|
writeRaster(scaledRaster,filename=outFile,format="HFA",overwrite=TRUE)
|
65
|
}
|
66
|
}
|
67
|
|
68
|
print("")
|
69
|
print("Processing h09v04 tiles")
|
70
|
print("")
|
71
|
|
72
|
# now reloop and process h09v04 tiles
|
73
|
# Create day averages for each day in the year
|
74
|
for(i in 1:366){
|
75
|
if (i < 9) {f <- "00"}
|
76
|
if ((i > 9) & (i <99 )) {f <- "0"}
|
77
|
if (i > 99) {f <- ""}
|
78
|
fpart <- paste(toString(f), toString(i), sep="")
|
79
|
|
80
|
# get a list of all files for that day of the year
|
81
|
print("")
|
82
|
print(sprintf("Reading day: %s", toString(fpart)))
|
83
|
myFiles <- grep(fpart,substr(inFiles, 14,16))
|
84
|
nFiles <- length(myFiles)
|
85
|
|
86
|
# for each file in the list
|
87
|
filestoStack <- list()
|
88
|
for (iCtr in 1 : nFiles)
|
89
|
{
|
90
|
mFile <- myFiles[iCtr]
|
91
|
if (substr(inFiles[mFile],18,23) == "h09v04") {
|
92
|
print(sprintf("..Reading file: %s",inFiles[mFile]))
|
93
|
inputRaster <- raster()
|
94
|
inputRaster <- raster(inFiles[mFile])
|
95
|
filestoStack <- c(filestoStack, inFiles[mFile])
|
96
|
}
|
97
|
}
|
98
|
# now calculate the mean cell values
|
99
|
myStack <- stack(filestoStack)
|
100
|
print(sprintf("....Averaging Daily Values"))
|
101
|
|
102
|
# Check if file exists
|
103
|
outFile <- paste("Day_",toString(i),"_Average_Scaled_h09v04",".img",sep="")
|
104
|
if (file.exists(outFile) == FALSE) {
|
105
|
meanRaster <- mean(myStack, na.rm=TRUE)
|
106
|
scaledRaster <- meanRaster * 0.0005 # apply scaling factor to raster
|
107
|
print(sprintf("......Writing file: %s",outFile))
|
108
|
writeRaster(scaledRaster,filename=outFile,format="HFA",overwrite=TRUE)
|
109
|
}
|
110
|
}
|
111
|
|
112
|
# finished
|
113
|
print("Finished")
|
114
|
|