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