1 |
147da66d
|
Adam M. Wilson
|
###################################################################################
|
2 |
|
|
### R code to aquire and process MOD06_L2 cloud data from the MODIS platform
|
3 |
|
|
|
4 |
|
|
|
5 |
|
|
## connect to server of choice
|
6 |
|
|
#system("ssh litoria")
|
7 |
|
|
#R
|
8 |
|
|
|
9 |
|
|
library(sp)
|
10 |
|
|
library(rgdal)
|
11 |
|
|
library(reshape)
|
12 |
|
|
library(ncdf4)
|
13 |
|
|
library(geosphere)
|
14 |
|
|
library(rgeos)
|
15 |
|
|
library(multicore)
|
16 |
|
|
library(raster)
|
17 |
|
|
library(lattice)
|
18 |
|
|
library(rgl)
|
19 |
|
|
library(hdf5)
|
20 |
|
|
library(spgrass6)
|
21 |
|
|
|
22 |
|
|
X11.options(type="Xlib")
|
23 |
|
|
ncores=20 #number of threads to use
|
24 |
|
|
|
25 |
|
|
setwd("/home/adamw/personal/projects/interp")
|
26 |
|
|
setwd("/home/adamw/acrobates/projects/interp")
|
27 |
|
|
|
28 |
|
|
roi=readOGR("data/regions/Test_sites/Oregon.shp","Oregon")
|
29 |
|
|
roi=spTransform(roi,CRS(" +proj=sinu +lon_0=0 +x_0=0 +y_0=0"))
|
30 |
|
|
|
31 |
|
|
### Downloading data from LAADSWEB
|
32 |
|
|
# subset by geographic area of interest
|
33 |
|
|
# subset: 40-47, -115--125
|
34 |
|
|
|
35 |
|
|
## download data from ftp site. Unfortunately the selection has to be selected using the website and orders downloaded via ftp.
|
36 |
|
|
|
37 |
|
|
system("wget -r --retr-symlinks ftp://ladsweb.nascom.nasa.gov/orders/500676499/ -P /home/adamw/acrobates/projects/interp/data/modis/MOD06_L2_hdf")
|
38 |
|
|
|
39 |
0c74b1da
|
Adam M. Wilson
|
## specify some working directories
|
40 |
147da66d
|
Adam M. Wilson
|
gdir="output/"
|
41 |
|
|
datadir="data/modis/MOD06_L2_hdf"
|
42 |
|
|
outdir="data/modis/MOD06_L2_hdf2"
|
43 |
0c74b1da
|
Adam M. Wilson
|
tifdir="/media/data/MOD06_L2_tif"
|
44 |
|
|
summarydatadir="data/modis/MOD06_climatologies"
|
45 |
|
|
|
46 |
08d8290d
|
Adam M. Wilson
|
|
47 |
147da66d
|
Adam M. Wilson
|
fs=data.frame(
|
48 |
|
|
path=list.files(datadir,full=T,recursive=T,pattern="hdf"),
|
49 |
|
|
file=basename(list.files(datadir,full=F,recursive=T,pattern="hdf")))
|
50 |
|
|
fs$date=as.Date(substr(fs$file,11,17),"%Y%j")
|
51 |
|
|
fs$month=format(fs$date,"%m")
|
52 |
|
|
fs$year=format(fs$date,"%Y")
|
53 |
|
|
fs$time=substr(fs$file,19,22)
|
54 |
|
|
fs$datetime=as.POSIXct(strptime(paste(substr(fs$file,11,17),substr(fs$file,19,22)), '%Y%j %H%M'))
|
55 |
979e2d4c
|
Adam M. Wilson
|
fs$dateid=format(fs$date,"%Y%m%d")
|
56 |
147da66d
|
Adam M. Wilson
|
fs$path=as.character(fs$path)
|
57 |
|
|
fs$file=as.character(fs$file)
|
58 |
|
|
|
59 |
|
|
## output ROI
|
60 |
|
|
#get bounding box of region in m
|
61 |
|
|
ge=SpatialPoints(data.frame(lon=c(-125,-115),lat=c(40,47)))
|
62 |
|
|
projection(ge)=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
|
63 |
|
|
ge2=spTransform(ge, CRS(" +proj=sinu +lon_0=0 +x_0=0 +y_0=0"))
|
64 |
|
|
|
65 |
|
|
|
66 |
|
|
## vars
|
67 |
|
|
vars=as.data.frame(matrix(c(
|
68 |
979e2d4c
|
Adam M. Wilson
|
"Cloud_Effective_Radius", "CER",
|
69 |
|
|
"Cloud_Effective_Radius_Uncertainty", "CERU",
|
70 |
|
|
"Cloud_Optical_Thickness", "COT",
|
71 |
|
|
"Cloud_Optical_Thickness_Uncertainty", "COTU",
|
72 |
|
|
"Cloud_Water_Path", "CWP",
|
73 |
|
|
"Cloud_Water_Path_Uncertainty", "CWPU",
|
74 |
|
|
"Cloud_Phase_Optical_Properties", "CPOP",
|
75 |
|
|
"Cloud_Multi_Layer_Flag", "CMLF",
|
76 |
|
|
"Cloud_Mask_1km", "CM1",
|
77 |
|
|
"Quality_Assurance_1km", "QA"),
|
78 |
0c74b1da
|
Adam M. Wilson
|
byrow=T,ncol=2,dimnames=list(1:10,c("variable","varid"))),stringsAsFactors=F)
|
79 |
147da66d
|
Adam M. Wilson
|
|
80 |
|
|
|
81 |
|
|
### Installation of hegtool
|
82 |
|
|
## needed 32-bit libraries and java for program to install correctly
|
83 |
|
|
|
84 |
|
|
# system(paste("hegtool -h ",fs$path[1],sep=""))
|
85 |
|
|
|
86 |
|
|
|
87 |
|
|
#### Function to generate hegtool parameter file for multi-band HDF-EOS file
|
88 |
0c74b1da
|
Adam M. Wilson
|
swath2grid=function(i=1,files,vars,outdir,upleft="47 -125",lowright="41 -115"){
|
89 |
147da66d
|
Adam M. Wilson
|
file=fs$path[i]
|
90 |
|
|
print(paste("Starting file",basename(file)))
|
91 |
|
|
outfile=paste(outdir,"/",fs$file[i],sep="")
|
92 |
|
|
# date=fs$date[1]
|
93 |
|
|
# origin=as.POSIXct("1970-01-01 00:00:00",tz="GMT")
|
94 |
|
|
### First write the parameter file (careful, heg is very finicky!)
|
95 |
|
|
hdr=paste("NUM_RUNS = ",length(vars),"|MULTI_BAND_HDFEOS:",length(vars),sep="")
|
96 |
|
|
grp=paste("
|
97 |
|
|
BEGIN
|
98 |
|
|
INPUT_FILENAME=",file,"
|
99 |
|
|
OBJECT_NAME=mod06
|
100 |
|
|
FIELD_NAME=",vars,"|
|
101 |
|
|
BAND_NUMBER = 1
|
102 |
|
|
OUTPUT_PIXEL_SIZE_X=1000
|
103 |
|
|
OUTPUT_PIXEL_SIZE_Y=1000
|
104 |
|
|
SPATIAL_SUBSET_UL_CORNER = ( ",upleft," )
|
105 |
|
|
SPATIAL_SUBSET_LR_CORNER = ( ",lowright," )
|
106 |
0c74b1da
|
Adam M. Wilson
|
#RESAMPLING_TYPE =",ifelse(grepl("Flag|Mask|Quality",vars),"NN","CUBIC"),"
|
107 |
147da66d
|
Adam M. Wilson
|
RESAMPLING_TYPE =NN
|
108 |
|
|
OUTPUT_PROJECTION_TYPE = SIN
|
109 |
|
|
ELLIPSOID_CODE = WGS84
|
110 |
979e2d4c
|
Adam M. Wilson
|
OUTPUT_PROJECTION_PARAMETERS = ( 6371007.181 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 86400.0 0.0 0.0 )
|
111 |
147da66d
|
Adam M. Wilson
|
OUTPUT_TYPE = HDFEOS
|
112 |
|
|
OUTPUT_FILENAME = ",outfile,"
|
113 |
|
|
END
|
114 |
|
|
|
115 |
|
|
",sep="")
|
116 |
|
|
## if any remnants from previous runs remain, delete them
|
117 |
|
|
if(length(list.files(tempdir(),pattern=basename(file)))>0)
|
118 |
|
|
file.remove(list.files(tempdir(),pattern=basename(file),full=T))
|
119 |
|
|
## write it to a file
|
120 |
|
|
cat(c(hdr,grp) , file=paste(tempdir(),"/",basename(file),"_MODparms.txt",sep=""))
|
121 |
|
|
## now run the swath2grid tool - must be run as root (argh!)!
|
122 |
|
|
## write the tiff file
|
123 |
|
|
log=system(paste("sudo MRTDATADIR=/usr/local/heg/data ",
|
124 |
|
|
"PGSHOME=/usr/local/heg/TOOLKIT_MTD PWD=/home/adamw /usr/local/heg/bin/swtif -p ",
|
125 |
|
|
paste(tempdir(),"/",basename(file),"_MODparms.txt -d",sep=""),sep=""),intern=T)
|
126 |
|
|
print(paste("Finished ", file))
|
127 |
|
|
}
|
128 |
08d8290d
|
Adam M. Wilson
|
|
129 |
147da66d
|
Adam M. Wilson
|
|
130 |
|
|
### update fs with completed files
|
131 |
|
|
fs$complete=fs$file%in%list.files(outdir,pattern="hdf$")
|
132 |
|
|
table(fs$complete)
|
133 |
|
|
|
134 |
|
|
#### Run the gridding procedure
|
135 |
|
|
|
136 |
|
|
system("sudo ls")
|
137 |
|
|
|
138 |
|
|
mclapply(which(!fs$complete),function(fi){
|
139 |
|
|
swath2grid(fi,vars=vars$variable,files=fs,
|
140 |
|
|
outdir=outdir,
|
141 |
|
|
upleft="47 -125",lowright="40 -115")},
|
142 |
|
|
mc.cores=24)
|
143 |
|
|
|
144 |
|
|
|
145 |
|
|
##############################################################
|
146 |
|
|
### Import to GRASS for processing
|
147 |
|
|
|
148 |
979e2d4c
|
Adam M. Wilson
|
#fs$grass=paste(fs$month,fs$year,fs$file,sep="_")
|
149 |
147da66d
|
Adam M. Wilson
|
td=readGDAL(paste("HDF4_EOS:EOS_GRID:\"",outdir,"/",fs$file[1],"\":mod06:Cloud_Mask_1km_0",sep=""))
|
150 |
0c74b1da
|
Adam M. Wilson
|
projection(td)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs +datum=WGS84 +ellps=WGS84 "
|
151 |
147da66d
|
Adam M. Wilson
|
|
152 |
|
|
## fucntion to convert binary to decimal to assist in identifying correct values
|
153 |
|
|
b2d=function(x) sum(x * 2^(rev(seq_along(x)) - 1)) #http://tolstoy.newcastle.edu.au/R/e2/help/07/02/10596.html
|
154 |
|
|
## for example:
|
155 |
|
|
b2d(c(T,T))
|
156 |
|
|
|
157 |
|
|
### create (or connect to) grass location
|
158 |
|
|
gisDbase="/media/data/grassdata"
|
159 |
|
|
gisLocation="oregon"
|
160 |
|
|
gisMapset="mod06"
|
161 |
08d8290d
|
Adam M. Wilson
|
## set Grass to overwrite
|
162 |
|
|
Sys.setenv(GRASS_OVERWRITE=1)
|
163 |
|
|
Sys.setenv(DEBUG=0)
|
164 |
147da66d
|
Adam M. Wilson
|
|
165 |
0c74b1da
|
Adam M. Wilson
|
## temporary objects to test function below
|
166 |
08d8290d
|
Adam M. Wilson
|
i=1
|
167 |
147da66d
|
Adam M. Wilson
|
file=paste(outdir,"/",fs$file[1],sep="")
|
168 |
979e2d4c
|
Adam M. Wilson
|
date=as.Date("2000-03-02")
|
169 |
|
|
|
170 |
0c74b1da
|
Adam M. Wilson
|
|
171 |
|
|
### Function to extract various SDSs from a single gridded HDF file and use QA data to throw out 'bad' observations
|
172 |
979e2d4c
|
Adam M. Wilson
|
loadcloud<-function(date,fs){
|
173 |
c33d3b68
|
Adam M. Wilson
|
### set up grass session
|
174 |
|
|
tf=paste(tempdir(),"/grass", Sys.getpid(),"/", sep="")
|
175 |
|
|
|
176 |
|
|
## set up tempfile for this PID
|
177 |
|
|
initGRASS(gisBase="/usr/lib/grass64",gisDbase=tf,SG=td,override=T,location="mod06",mapset="PERMANENT",home=tf,pid=Sys.getpid())
|
178 |
|
|
system(paste("g.proj -c proj4=\"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +datum=WGS84 +units=m +no_defs\"",sep=""))
|
179 |
|
|
|
180 |
|
|
#system("g.mapset PERMANENT")
|
181 |
|
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",outdir,"/",fs$file[1],"\":mod06:Cloud_Mask_1km_0",sep=""),
|
182 |
|
|
output="modisgrid",flags=c("quiet","overwrite","o"))
|
183 |
|
|
system("g.region rast=modisgrid save=roi --overwrite")
|
184 |
|
|
system("g.region roi")
|
185 |
|
|
system("g.region -p")
|
186 |
|
|
# getLocationProj()
|
187 |
|
|
|
188 |
|
|
|
189 |
|
|
## Identify which files to process
|
190 |
979e2d4c
|
Adam M. Wilson
|
tfs=fs$file[fs$date==date]
|
191 |
|
|
nfs=length(tfs)
|
192 |
c33d3b68
|
Adam M. Wilson
|
|
193 |
|
|
### print some summary info
|
194 |
08d8290d
|
Adam M. Wilson
|
print(date)
|
195 |
979e2d4c
|
Adam M. Wilson
|
## loop through scenes and process QA flags
|
196 |
|
|
for(i in 1:nfs){
|
197 |
08d8290d
|
Adam M. Wilson
|
file=paste(outdir,"/",tfs[i],sep="")
|
198 |
|
|
## Cloud Mask
|
199 |
|
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Mask_1km_0",sep=""),
|
200 |
|
|
output=paste("CM1_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
201 |
979e2d4c
|
Adam M. Wilson
|
## extract cloudy and 'confidently clear' pixels
|
202 |
|
|
system(paste("r.mapcalc <<EOF
|
203 |
|
|
CM_cloud_",i," = ((CM1_",i," / 2^0) % 2) == 1 && ((CM1_",i," / 2^1) % 2^2) == 0
|
204 |
|
|
CM_clear_",i," = ((CM1_",i," / 2^0) % 2) == 1 && ((CM1_",i," / 2^1) % 2^2) == 3
|
205 |
|
|
EOF",sep=""))
|
206 |
|
|
|
207 |
|
|
## QA
|
208 |
08d8290d
|
Adam M. Wilson
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Quality_Assurance_1km_0",sep=""),
|
209 |
|
|
output=paste("QA_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
210 |
147da66d
|
Adam M. Wilson
|
## QA_CER
|
211 |
|
|
system(paste("r.mapcalc <<EOF
|
212 |
979e2d4c
|
Adam M. Wilson
|
QA_COT_",i,"= ((QA_",i," / 2^0) % 2^1 )==1
|
213 |
|
|
QA_COT2_",i,"= ((QA_",i," / 2^1) % 2^2 )==3
|
214 |
|
|
QA_COT3_",i,"= ((QA_",i," / 2^3) % 2^2 )==0
|
215 |
|
|
QA_CER_",i,"= ((QA_",i," / 2^5) % 2^1 )==1
|
216 |
|
|
QA_CER2_",i,"= ((QA_",i," / 2^6) % 2^2 )==3
|
217 |
|
|
EOF",sep=""))
|
218 |
c33d3b68
|
Adam M. Wilson
|
# QA_CWP_",i,"= ((QA_",i," / 2^8) % 2^1 )==1
|
219 |
|
|
# QA_CWP2_",i,"= ((QA_",i," / 2^9) % 2^2 )==3
|
220 |
147da66d
|
Adam M. Wilson
|
|
221 |
|
|
## Optical Thickness
|
222 |
|
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Optical_Thickness",sep=""),
|
223 |
979e2d4c
|
Adam M. Wilson
|
output=paste("COT_",i,sep=""),
|
224 |
147da66d
|
Adam M. Wilson
|
title="cloud_effective_radius",
|
225 |
|
|
flags=c("overwrite","o")) ; print("")
|
226 |
979e2d4c
|
Adam M. Wilson
|
execGRASS("r.null",map=paste("COT_",i,sep=""),setnull="-9999")
|
227 |
147da66d
|
Adam M. Wilson
|
## keep only positive COT values where quality is 'useful' and 'very good' & scale to real units
|
228 |
979e2d4c
|
Adam M. Wilson
|
system(paste("r.mapcalc \"COT_",i,"=if(QA_COT_",i,"&&QA_COT2_",i,"&&QA_COT3_",i,"&&COT_",i,">=0,COT_",i,"*0.009999999776482582,null())\"",sep=""))
|
229 |
147da66d
|
Adam M. Wilson
|
## set COT to 0 in clear-sky pixels
|
230 |
979e2d4c
|
Adam M. Wilson
|
system(paste("r.mapcalc \"COT2_",i,"=if(CM_clear_",i,"==0,COT_",i,",0)\"",sep=""))
|
231 |
147da66d
|
Adam M. Wilson
|
|
232 |
979e2d4c
|
Adam M. Wilson
|
## Effective radius ##
|
233 |
147da66d
|
Adam M. Wilson
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Effective_Radius",sep=""),
|
234 |
979e2d4c
|
Adam M. Wilson
|
output=paste("CER_",i,sep=""),
|
235 |
147da66d
|
Adam M. Wilson
|
title="cloud_effective_radius",
|
236 |
|
|
flags=c("overwrite","o")) ; print("")
|
237 |
979e2d4c
|
Adam M. Wilson
|
execGRASS("r.null",map=paste("CER_",i,sep=""),setnull="-9999")
|
238 |
147da66d
|
Adam M. Wilson
|
## keep only positive CER values where quality is 'useful' and 'very good' & scale to real units
|
239 |
979e2d4c
|
Adam M. Wilson
|
system(paste("r.mapcalc \"CER_",i,"=if(QA_CER_",i,"&&QA_CER2_",i,"&&CER_",i,">=0,CER_",i,"*0.009999999776482582,null())\"",sep=""))
|
240 |
|
|
## set CER to 0 in clear-sky pixels
|
241 |
|
|
system(paste("r.mapcalc \"CER2_",i,"=if(CM_clear_",i,"==0,CER_",i,",0)\"",sep=""))
|
242 |
147da66d
|
Adam M. Wilson
|
|
243 |
|
|
## Cloud Water Path
|
244 |
c33d3b68
|
Adam M. Wilson
|
# execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Water_Path",sep=""),
|
245 |
|
|
# output=paste("CWP_",i,sep=""),title="cloud_water_path",
|
246 |
|
|
# flags=c("overwrite","o")) ; print("")
|
247 |
|
|
# execGRASS("r.null",map=paste("CWP_",i,sep=""),setnull="-9999")
|
248 |
0c74b1da
|
Adam M. Wilson
|
## keep only positive CER values where quality is 'useful' and 'very good' & scale to real units
|
249 |
c33d3b68
|
Adam M. Wilson
|
# system(paste("r.mapcalc \"CWP_",i,"=if(QA_CWP_",i,"&&QA_CWP2_",i,"&&CWP_",i,">=0,CWP_",i,"*0.009999999776482582,null())\"",sep=""))
|
250 |
0c74b1da
|
Adam M. Wilson
|
## set CER to 0 in clear-sky pixels
|
251 |
c33d3b68
|
Adam M. Wilson
|
# system(paste("r.mapcalc \"CWP2_",i,"=if(CM_clear_",i,"==0,CWP_",i,",0)\"",sep=""))
|
252 |
0c74b1da
|
Adam M. Wilson
|
|
253 |
|
|
|
254 |
979e2d4c
|
Adam M. Wilson
|
} #end loop through sub daily files
|
255 |
|
|
|
256 |
0c74b1da
|
Adam M. Wilson
|
#### Now generate daily averages (or maximum in case of cloud flag)
|
257 |
|
|
|
258 |
979e2d4c
|
Adam M. Wilson
|
system(paste("r.mapcalc <<EOF
|
259 |
|
|
COT_denom=",paste("!isnull(COT2_",1:nfs,")",sep="",collapse="+"),"
|
260 |
|
|
COT_numer=",paste("if(isnull(COT2_",1:nfs,"),0,COT2_",1:nfs,")",sep="",collapse="+"),"
|
261 |
|
|
COT_daily=COT_numer/COT_denom
|
262 |
|
|
CER_denom=",paste("!isnull(CER2_",1:nfs,")",sep="",collapse="+"),"
|
263 |
|
|
CER_numer=",paste("if(isnull(CER2_",1:nfs,"),0,CER2_",1:nfs,")",sep="",collapse="+"),"
|
264 |
|
|
CER_daily=CER_numer/CER_denom
|
265 |
0c74b1da
|
Adam M. Wilson
|
CLD_daily=max(",paste("if(isnull(CM_cloud_",1:nfs,"),0,CM_cloud_",1:nfs,")",sep="",collapse=","),")
|
266 |
979e2d4c
|
Adam M. Wilson
|
EOF",sep=""))
|
267 |
|
|
|
268 |
|
|
#### Write the file to a geotiff
|
269 |
c33d3b68
|
Adam M. Wilson
|
execGRASS("r.out.gdal",input="CER_daily",output=paste(tifdir,"/CER_",format(date,"%Y%m%d"),".tif",sep=""),nodata=-999,flags=c("quiet"))
|
270 |
|
|
execGRASS("r.out.gdal",input="COT_daily",output=paste(tifdir,"/COT_",format(date,"%Y%m%d"),".tif",sep=""),nodata=-999,flags=c("quiet"))
|
271 |
|
|
execGRASS("r.out.gdal",input="CLD_daily",output=paste(tifdir,"/CLD_",format(date,"%Y%m%d"),".tif",sep=""),nodata=99,flags=c("quiet"))
|
272 |
08d8290d
|
Adam M. Wilson
|
|
273 |
|
|
### delete the temporary files
|
274 |
|
|
unlink_.gislock()
|
275 |
|
|
system("/usr/lib/grass64/etc/clean_temp")
|
276 |
c33d3b68
|
Adam M. Wilson
|
system(paste("rm -R ",tf,sep=""))
|
277 |
|
|
### print update
|
278 |
|
|
print(paste(" ################################################################### Finished ",date,"
|
279 |
|
|
################################################################"))
|
280 |
979e2d4c
|
Adam M. Wilson
|
}
|
281 |
|
|
|
282 |
|
|
|
283 |
|
|
###########################################
|
284 |
|
|
### Now run it
|
285 |
|
|
|
286 |
0c74b1da
|
Adam M. Wilson
|
tdates=sort(unique(fs$date))
|
287 |
c33d3b68
|
Adam M. Wilson
|
done=tdates%in%as.Date(substr(list.files(tifdir),5,12),"%Y%m%d")
|
288 |
08d8290d
|
Adam M. Wilson
|
table(done)
|
289 |
|
|
tdates=tdates[!done]
|
290 |
979e2d4c
|
Adam M. Wilson
|
|
291 |
c33d3b68
|
Adam M. Wilson
|
mclapply(tdates,function(date) loadcloud(date,fs=fs))
|
292 |
147da66d
|
Adam M. Wilson
|
|
293 |
0c74b1da
|
Adam M. Wilson
|
|
294 |
147da66d
|
Adam M. Wilson
|
|
295 |
08d8290d
|
Adam M. Wilson
|
#######################################################################################33
|
296 |
|
|
### Produce the monthly averages
|
297 |
147da66d
|
Adam M. Wilson
|
|
298 |
08d8290d
|
Adam M. Wilson
|
## get list of daily files
|
299 |
|
|
fs2=data.frame(
|
300 |
0c74b1da
|
Adam M. Wilson
|
path=list.files(tifdir,full=T,recursive=T,pattern="tif$"),
|
301 |
|
|
file=basename(list.files(tifdir,full=F,recursive=T,pattern="tif$")))
|
302 |
|
|
fs2$type=substr(fs2$file,1,3)
|
303 |
|
|
fs2$date=as.Date(substr(fs2$file,5,12),"%Y%m%d")
|
304 |
|
|
fs2$month=format(fs2$date,"%m")
|
305 |
|
|
fs2$year=format(fs2$date,"%Y")
|
306 |
|
|
fs2$path=as.character(fs2$path)
|
307 |
|
|
fs2$file=as.character(fs2$file)
|
308 |
|
|
|
309 |
|
|
|
310 |
|
|
# Define type/month products
|
311 |
|
|
vs=expand.grid(type=unique(fs2$type),month=c("01","02","03","04","05","06","07","08","09","10","11","12"))
|
312 |
|
|
|
313 |
|
|
## identify which have been completed
|
314 |
|
|
#done=
|
315 |
|
|
# do.call(rbind,strsplit(list.files(summarydatadir),"_|[.]"))[,3]
|
316 |
|
|
#table(done)
|
317 |
|
|
#tdates=tdates[!done]
|
318 |
|
|
|
319 |
|
|
|
320 |
|
|
## process the summaries using the raster package
|
321 |
|
|
mclapply(1:nrow(vs),function(i){
|
322 |
|
|
print(paste("Starting ",vs$type[i]," for month ",vs$month[i]))
|
323 |
|
|
td=stack(fs2$path[which(fs2$month==vs$month[i]&fs2$type==vs$type[i])])
|
324 |
c33d3b68
|
Adam M. Wilson
|
print(paste("Processing Metric ",vs$type[i]," for month ",vs$month[i]))
|
325 |
0c74b1da
|
Adam M. Wilson
|
calc(td,mean,na.rm=T,
|
326 |
|
|
filename=paste(summarydatadir,"/",vs$type[i],"_mean_",vs$month[i],".tif",sep=""),
|
327 |
|
|
format="GTiff")
|
328 |
|
|
calc(td,sd,na.rm=T,
|
329 |
|
|
filename=paste(summarydatadir,"/",vs$type[i],"_sd_",vs$month[i],".tif",sep=""),
|
330 |
|
|
format="GTiff")
|
331 |
|
|
print(paste("Processing missing data for ",vs$type[i]," for month ",vs$month[i]))
|
332 |
|
|
calc(td,function(i)
|
333 |
|
|
sum(!is.na(i)),filename=paste(summarydatadir,"/",vs$type[i],"_count_",vs$month[i],".tif",sep=""),
|
334 |
|
|
format="GTiff")
|
335 |
|
|
calc(td,function(i) sum(ifelse(i==0,0,1)),
|
336 |
|
|
filename=paste(summarydatadir,"/",vs$type[i],"_clear_",vs$month[i],".tif",sep=""),format="GTiff")
|
337 |
147da66d
|
Adam M. Wilson
|
gc()
|
338 |
|
|
}
|
339 |
0c74b1da
|
Adam M. Wilson
|
)
|
340 |
147da66d
|
Adam M. Wilson
|
|