Revision 31bee7b1
Added by Adam Wilson almost 11 years ago
| climate/research/cloud/MOD09/2_bias.R | ||
|---|---|---|
| 4 | 4 |
library(doMC) |
| 5 | 5 |
library(foreach) |
| 6 | 6 |
library(RcppOctave) |
| 7 |
registerDoMC(12) |
|
| 7 |
library(rgdal) |
|
| 8 |
registerDoMC(7) |
|
| 8 | 9 |
|
| 9 | 10 |
|
| 10 |
## start raster cluster |
|
| 11 |
#beginCluster(5) |
|
| 12 |
|
|
| 11 |
# final output will be written to data directory here: |
|
| 13 | 12 |
setwd("~/acrobates/adamw/projects/cloud")
|
| 14 | 13 |
|
| 14 |
# temporary files will be written here: |
|
| 15 | 15 |
datadir="/mnt/data2/projects/cloud/" |
| 16 | 16 |
|
| 17 |
|
|
| 18 |
## Specify path to VSNR souce code and add it to RcppOctave path |
|
| 17 | 19 |
mpath="/home/adamw/acrobates/adamw/projects/environmental-layers/climate/research/cloud/MOD09/vsnr/" |
| 18 | 20 |
.O$addpath(mpath) |
| 19 | 21 |
|
| ... | ... | |
| 39 | 41 |
} |
| 40 | 42 |
|
| 41 | 43 |
|
| 42 |
vsnr=function(d,gabor,alpha=10,p=2,epsilon=0,prec=5e-3,maxit=100,C1=1){
|
|
| 43 |
# d=d*1.0 |
|
| 44 |
vsnr=function(d,gabor,alpha=2,p=2,epsilon=0,prec=5e-3,maxit=100,C1=1,full=F){
|
|
| 45 |
## VSNR can't run with any missing values, set them to zero here then switch them back to NA later |
|
| 46 |
d2=as.matrix(d) |
|
| 47 |
d2[is.na(d2)]=0 |
|
| 44 | 48 |
## Process with VSNR |
| 45 |
dt=.O$VSNR(as.matrix(d),epsilon,p,as.matrix(gabor),alpha,maxit+1,prec,1,argout = c("u", "Gap", "Primal","Dual","EstP","EstD"));
|
|
| 49 |
dt=.CallOctave("VSNR",d2,epsilon,p,as.matrix(gabor),alpha,maxit+1,prec,C1,argout ="u",verbose=T);
|
|
| 46 | 50 |
## Create spatial objects from VSNR output |
| 47 |
bias=d |
|
| 48 |
values(bias)=-as.numeric(t(convolve(dt$EstP,as.matrix(gabor)))) |
|
| 49 |
#bias[abs(bias)<5]=0 # set all bias<0.5 to zero to avoit minor adjustment and added striping in data |
|
| 50 |
#bias=bias*10 |
|
| 51 |
# dc=d-bias # Make 'corrected' version of data |
|
| 52 |
# NAvalue(bias)=0 #set bias==0 to NA to facilate plotting, areas that are NA were not adjusted |
|
| 53 |
# res=stack(list(dc=dc,d=d,bias=bias)) #,EstP=EstP,EstD1=EstD1,EstD2=EstD2 |
|
| 54 |
res=stack(list(bias=bias)) #,EstP=EstP,EstD1=EstD1,EstD2=EstD2 |
|
| 55 |
# rm(dt,dc) |
|
| 56 |
rm(dt) |
|
| 57 |
return(res) |
|
| 51 |
dc=d |
|
| 52 |
values(dc)=as.numeric(t(dt$u)) # Make 'corrected' version of data |
|
| 53 |
## Set NA values in original data back to NA |
|
| 54 |
dc[is.na(d)]=NA |
|
| 55 |
return(dc) |
|
| 58 | 56 |
} |
| 59 | 57 |
|
| 58 |
rmr=function(x){
|
|
| 59 |
## function to truly delete raster and temporary files associated with them |
|
| 60 |
if(class(x)=="RasterLayer"&grepl("^/tmp",x@file@name)&fromDisk(x)==T){
|
|
| 61 |
file.remove(x@file@name,sub("grd","gri",x@file@name))
|
|
| 62 |
rm(x) |
|
| 63 |
} |
|
| 64 |
} |
|
| 60 | 65 |
|
| 66 |
|
|
| 61 | 67 |
###################################### |
| 62 | 68 |
## Run the correction functions |
| 63 | 69 |
### Subset equitorial region to correct orbital banding |
| ... | ... | |
| 75 | 81 |
|
| 76 | 82 |
df2=list.files(paste(datadir,"/mcd09tif",sep=""),full=T,pattern="[0-9][.]tif$") |
| 77 | 83 |
i=1 |
| 84 |
ti=7 |
|
| 78 | 85 |
|
| 79 | 86 |
### Build the tiles |
| 80 |
exts=extf(xmin=-180,xmax=180,ymin=-30,ymax=30,size=10,overlap=0) |
|
| 81 |
#exts=extf(xmin=-10,xmax=20,ymin=0,ymax=30,size=10,overlap=0) |
|
| 82 |
#exts=extf(xmin=-60,xmax=60,ymin=-30,ymax=30,size=10,overlap=0) |
|
| 87 |
exts=extf(xmin=-180,xmax=180,ymin=-30,ymax=30,size=60,overlap=0) |
|
| 88 |
|
|
| 89 |
## add an extra tile to account for regions of reduced data availability for each sensor |
|
| 90 |
modexts=c(xmin=130,xmax=180,ymin=-50,ymax=0) |
|
| 91 |
mydexts=c(xmin=-170,xmax=-140,ymin=-30,ymax=55) |
|
| 83 | 92 |
|
| 84 | 93 |
|
| 85 | 94 |
## loop over sensor-months to create full grid of corrected values |
| 86 |
for( i in 1:length(df2)){
|
|
| 87 |
### Process the tiles |
|
| 88 |
foreach(ti=1:nrow(exts)) %dopar% {
|
|
| 89 |
textent=extent(exts$xmin[ti],exts$xmax[ti],exts$ymin[ti],exts$ymax[ti]) |
|
| 90 |
## extract the tile |
|
| 95 |
for( i in 1:nrow(df2))){
|
|
| 91 | 96 |
file=df2[i] |
| 92 |
outfile=paste("tmp/",sub(".tif","",basename(file)),"_",sprintf("%03d",ti),".tif",sep="")
|
|
| 93 |
writeLines(paste("Starting: ",outfile," tile:",ti," ( out of ",nrow(exts),")"))
|
|
| 94 |
d=crop(raster(file),textent) |
|
| 95 |
## acount for scale of data is 10000*CF |
|
| 96 |
d=d*.01 |
|
| 97 |
## skip null tiles |
|
| 98 |
if(is.null(d@data@values)) return(NULL) |
|
| 99 |
## make the gabor kernel |
|
| 100 |
psi=fgabor(d,theta=-15,x=400,y=4) #3 |
|
| 101 |
## run the correction |
|
| 102 |
res=vsnr(d,gabor=psi,alpha=2,p=2,epsilon=1,prec=5e-6,maxit=50,C=1) |
|
| 103 |
## write the file |
|
| 104 |
if(!is.null(outfile)) writeRaster(res,file=outfile,overwrite=T,datatype='INT2S',options=c("COMPRESS=LZW", "PREDICTOR=2"))#,NAflag=-127)
|
|
| 105 |
print(paste("Finished ",outfile))
|
|
| 106 |
} |
|
| 97 |
outfile=paste(datadir,"/mcd09ctif/",basename(file),sep="") |
|
| 98 |
# outfile2=paste("data/mcd09tif/",basename(file),sep="")
|
|
| 99 |
# outfile2b=paste("data/mcd09tif/",sub("[.]tif","_sd.tif",basename(file)),sep="")
|
|
| 100 |
|
|
| 101 |
## set sensor-specific parameters |
|
| 102 |
## add extra region for correction depending on which sensor is being processed |
|
| 103 |
## set angle of orbital artefacts to be corrected |
|
| 104 |
sensor=ifelse(grepl("MOD",file),"MOD","MYD")
|
|
| 105 |
if(sensor=="MOD") {
|
|
| 106 |
exts=rbind(exts,modexts) |
|
| 107 |
scanangle=-15 |
|
| 108 |
} |
|
| 109 |
if(sensor=="MYD") {
|
|
| 110 |
exts=rbind(exts,mydexts) |
|
| 111 |
scanangle=15 |
|
| 112 |
} |
|
| 113 |
|
|
| 114 |
## Process the tiles |
|
| 115 |
foreach(ti=1:nrow(exts)) %dopar% {
|
|
| 116 |
textent=extent(exts$xmin[ti],exts$xmax[ti],exts$ymin[ti],exts$ymax[ti]) |
|
| 117 |
## extract the tile |
|
| 118 |
toutfile=paste("tmp/", sub(".tif","",basename(file)),"_",sprintf("%03d",ti),".tif",sep="")
|
|
| 119 |
writeLines(paste("Starting: ",toutfile," tile:",ti," ( out of ",nrow(exts),")"))
|
|
| 120 |
d=crop(raster(file),textent) |
|
| 121 |
## acount for scale of data is 10000*CF |
|
| 122 |
d=d*.01 |
|
| 123 |
## skip null tiles - will only have this if tiles are quite small (<10 degrees) |
|
| 124 |
if(is.null(d@data@values)) return(NULL) |
|
| 125 |
## make the gabor kernel |
|
| 126 |
## this specifies the 'shape' of the noise we are trying to remove |
|
| 127 |
psi=fgabor(d,theta=scanangle,x=400,y=4) #3 |
|
| 128 |
# psi=stack(lapply(scanangle,function(a) fgabor(d,theta=a,x=400,y=4))) #3 |
|
| 129 |
|
|
| 130 |
## run the correction function. |
|
| 131 |
res=vsnr(d,gabor=psi,alpha=2,p=2,epsilon=1,prec=5e-6,maxit=50,C=1,full=F) |
|
| 132 |
## write the file |
|
| 133 |
if(!is.null(outfile)) writeRaster(res*100,file=toutfile,overwrite=T,datatype='INT2S',options=c("COMPRESS=LZW", "PREDICTOR=2"),NAvalue=-32768)
|
|
| 134 |
## remove temporary files |
|
| 135 |
rmr(d);rmr(psi);rmr(res) |
|
| 136 |
print(paste("Finished Temporary File: ",toutfile))
|
|
| 137 |
} |
|
| 138 |
## create VRT of first band of the full image |
|
| 139 |
fvrt=sub("[.]tif","_cf.vrt",file)
|
|
| 140 |
system(paste("gdalbuildvrt -b 1 ",fvrt," ",file))
|
|
| 141 |
## mosaic the tiles with the original data (keeping the new data when available) |
|
| 142 |
tfiles=paste(c(fvrt,list.files("tmp",pattern=paste(sub("[.]tif","",basename(outfile)),"_[0-9]*[.]tif",sep=""),full=T)),collapse=" ")
|
|
| 143 |
system(paste("gdal_merge.py -init -32768 -n -32768 -co COMPRESS=LZW -co PREDICTOR=2 -co BIGTIFF=yes -o ",outfile," ",tfiles,sep=""))
|
|
| 107 | 144 |
} |
| 108 | 145 |
|
| 109 |
## mosaic the tiles |
|
| 110 |
system("ls tmp/test_[0-9]*[.]tif")
|
|
| 111 |
system("rm tmp/test2.tif; gdal_merge.py -co COMPRESS=LZW -co PREDICTOR=2 -co BIGTIFF=yes 'tmp/test_[0-9]*[.]tif' -o tmp/test2.tif")
|
|
| 112 |
#system("rm tmp/test2.tif; gdalwarp -multi -r average -dstnodata 255 -ot Byte -co COMPRESS=LZW -co PREDICTOR=2 `ls tmp/test_[0-9]*[.]tif` tmp/test2.tif")
|
|
| 113 |
|
|
| 114 |
res=brick("tmp/test2.tif")
|
|
| 115 |
names(res)=c("dc","d","bias")
|
|
| 116 |
NAvalue(res)=-32768 |
|
| 117 |
|
|
| 118 |
tplot=F |
|
| 119 |
if(tplot){
|
|
| 120 |
gcol=colorRampPalette(c("blue","yellow","red"))
|
|
| 121 |
gcol=colorRampPalette(c("black","white"))
|
|
| 122 |
levelplot(res[["bias"]],col.regions=gcol(100),cuts=99,margin=F,maxpixels=1e6) |
|
| 123 |
levelplot(stack(list(Original=res[["d"]],Corrected=res[["dc"]])),col.regions=gcol(100),cuts=99,maxpixels=1e6) |
|
| 124 |
levelplot(stack(list(Original=res[["d"]],Corrected=res[["dc"]])),col.regions=gcol(100),cuts=99,maxpixels=1e6,ylim=c(10,13),xlim=c(-7,-1)) |
|
| 125 |
} |
|
| 146 |
|
|
| 147 |
################################################################################ |
|
| 148 |
### calculate monthly means of terra and aqua |
|
| 149 |
|
|
| 150 |
# Create combined (MOD+MYD) uncorrected mean CF |
|
| 151 |
foreach(i=1:12) %dopar% {
|
|
| 152 |
## get files |
|
| 153 |
f=list.files(paste(datadir,"/mcd09ctif",sep=""),pattern=paste(".*[O|Y].*_",sprintf("%02d",i),"[.]tif$",sep=""),full=T)
|
|
| 154 |
## Define output and check if it already exists |
|
| 155 |
tmcd=paste(datadir,"/mcd09ctif/MCD09_",sprintf("%02d", i),".tif",sep="")
|
|
| 156 |
## check if output already exists |
|
| 157 |
ops=paste("-t_srs 'EPSG:4326' -multi -srcnodata -32768 -dstnodata -32768 -r bilinear -te -180 -90 180 90 -tr 0.008333333333333 -0.008333333333333",
|
|
| 158 |
"-co BIGTIFF=YES --config GDAL_CACHEMAX 500 -wm 500 -wo NUM_THREADS:10 -wo SOURCE_EXTRA=5") |
|
| 159 |
system(paste("gdalwarp -overwrite -r average -co COMPRESS=LZW -co ZLEVEL=9 ",ops," ",paste(f,collapse=" ")," ",tmcd))
|
|
| 160 |
## update metadata |
|
| 161 |
tags=c(paste("TIFFTAG_IMAGEDESCRIPTION='Monthly Cloud Frequency for 2000-2013 extracted from C5 MODIS MOD09GA and MYD09GA PGE11 internal cloud mask algorithm (embedded in state_1km bit 10).",
|
|
| 162 |
"The daily cloud mask time series were summarized to mean cloud frequency (CF) by calculating the proportion of cloudy days. ", |
|
| 163 |
"Band Descriptions: 1) Mean Monthly Cloud Frequency 2) Four Times the Standard Deviation of Mean Monthly Cloud Frequency"'"), |
|
| 164 |
"TIFFTAG_DOCUMENTNAME='Collection 5 MCD09 Cloud Frequency'", |
|
| 165 |
paste("TIFFTAG_DATETIME='2013",sprintf("%02d", i),"15'",sep=""),
|
|
| 166 |
"TIFFTAG_ARTIST='Adam M. Wilson (adam.wilson@yale.edu)'") |
|
| 167 |
system(paste("/usr/local/src/gdal-1.10.0/swig/python/scripts/gdal_edit.py ",tmcd," ",paste("-mo ",tags,sep="",collapse=" "),sep=""))
|
|
| 168 |
writeLines(paste("Finished month",i))
|
|
| 169 |
} |
|
| 170 |
|
|
| 126 | 171 |
|
| 127 | 172 |
|
| climate/research/cloud/MOD09/3_combine.R | ||
|---|---|---|
| 1 |
################################################################################ |
|
| 2 |
### calculate monthly means of terra and aqua |
|
| 3 |
|
|
| 4 |
datadir="/mnt/data2/projects/cloud/" |
|
| 5 |
|
|
| 6 |
# Create combined (MOD+MYD) corrected mean CF |
|
| 7 |
foreach(i=1:12) %dopar% {
|
|
| 8 |
|
|
| 9 |
f=list.files(paste(datadir,"/mcd09ctif",sep=""),pattern=paste(".*[O|Y].*_",sprintf("%02d",i),"[.]tif$",sep=""),full=T)
|
|
| 10 |
## Define output and check if it already exists |
|
| 11 |
tmcd=paste(datadir,"/mcd09ctif/MCD09_",sprintf("%02d", i),".tif",sep="")
|
|
| 12 |
## check if output already exists |
|
| 13 |
ops=paste("-t_srs 'EPSG:4326' -multi -srcnodata -32768 -dstnodata -32768 -r bilinear -te -180 -90 180 90 -tr 0.008333333333333 -0.008333333333333",
|
|
| 14 |
"-co BIGTIFF=YES --config GDAL_CACHEMAX 500 -wm 500 -wo NUM_THREADS:10 -wo SOURCE_EXTRA=5") |
|
| 15 |
system(paste("gdalwarp -overwrite -r average -co COMPRESS=LZW -co ZLEVEL=9 ",ops," ",paste(f,collapse=" ")," ",tmcd))
|
|
| 16 |
## update metadata |
|
| 17 |
tags=c(paste("TIFFTAG_IMAGEDESCRIPTION='Monthly Cloud Frequency for 2000-2013 extracted from C5 MODIS MOD09GA and MYD09GA PGE11 internal cloud mask algorithm (embedded in state_1km bit 10).",
|
|
| 18 |
"The daily cloud mask time series were summarized to mean cloud frequency (CF) by calculating the proportion of cloudy days. ", |
|
| 19 |
"Band Descriptions: 1) Mean Monthly Cloud Frequency 2) Four Times the Standard Deviation of Mean Monthly Cloud Frequency", |
|
| 20 |
" 3) Mean number of daily observations for each pixel 4) Proportion of days with at least one observation '"), |
|
| 21 |
"TIFFTAG_DOCUMENTNAME='Collection 5 MCD09 Cloud Frequency'", |
|
| 22 |
paste("TIFFTAG_DATETIME='2013",sprintf("%02d", i),"15'",sep=""),
|
|
| 23 |
"TIFFTAG_ARTIST='Adam M. Wilson (adam.wilson@yale.edu)'") |
|
| 24 |
system(paste("/usr/local/src/gdal-1.10.0/swig/python/scripts/gdal_edit.py ",tmcd," ",paste("-mo ",tags,sep="",collapse=" "),sep=""))
|
|
| 25 |
writeLines(paste("Finished month",i))
|
|
| 26 |
} |
|
| 27 |
|
|
| 28 |
|
|
| 29 |
|
|
| 30 |
################################################################################# |
|
| 31 |
###### convert to 8-bit compressed file, add colors and other details |
|
| 32 |
for( i in 1:12){
|
|
| 33 |
f2=list.files(paste(datadir,"/mcd09ctif",sep=""),pattern=paste(".*MCD09_",sprintf("%02d",i),"[.]tif$",sep=""),full=T)
|
|
| 34 |
|
|
| 35 |
outfile=paste(datadir,"/mcd09ctif/",basename(file),sep="") |
|
| 36 |
outfile2=paste("data/mcd09tif/MCD09_",sprintf("%02d",i),".tif",sep="")
|
|
| 37 |
outfile2b=paste("data/mcd09tif/MCD09_",sprintf("%02d",i),"_sd.tif",sep="")
|
|
| 38 |
|
|
| 39 |
## convert to 8-bit compressed file |
|
| 40 |
ops="-stats -scale 0 10000 0 100 -ot Byte -a_nodata 255 -co COMPRESS=LZW -co PREDICTOR=2 -co BIGTIFF=yes" |
|
| 41 |
tags=c(paste("TIFFTAG_IMAGEDESCRIPTION='Monthly Cloud Frequency for 2000-2013 extracted from C5 MODIS M*D09GA PGE11 internal cloud mask algorithm (embedded in state_1km bit 10).",
|
|
| 42 |
"The daily cloud mask time series were summarized to mean cloud frequency (CF) by calculating the proportion of cloudy days. ", |
|
| 43 |
"Band Descriptions: 1) Mean Monthly Cloud Frequency'"), |
|
| 44 |
"TIFFTAG_DOCUMENTNAME='Collection 5 Cloud Frequency'", |
|
| 45 |
paste("TIFFTAG_DATETIME='2014'",sep=""),
|
|
| 46 |
"TIFFTAG_ARTIST='Adam M. Wilson (adam.wilson@yale.edu)'") |
|
| 47 |
system(paste("gdal_translate ",ops," ",paste("-mo ",tags,sep="",collapse=" ")," ",outfile," ",outfile2))
|
|
| 48 |
## Convert SD to 8-bit |
|
| 49 |
system(paste("gdal_translate -b 2 ",ops," ",paste("-mo ",tags,sep="",collapse=" ")," ",file," ",outfile2b))
|
|
| 50 |
} |
|
| 51 |
|
|
| 52 |
|
|
| 53 |
tplot=F |
|
| 54 |
if(tplot){
|
|
| 55 |
gcol=colorRampPalette(c("blue","yellow","red"))
|
|
| 56 |
gcol=colorRampPalette(c("black","white"))
|
|
| 57 |
levelplot(res[["bias"]],col.regions=gcol(100),cuts=99,margin=F,maxpixels=1e6) |
|
| 58 |
levelplot(stack(list(Original=res[["d"]],Corrected=res[["dc"]])),col.regions=gcol(100),cuts=99,maxpixels=1e6) |
|
| 59 |
levelplot(stack(list(Original=res[["d"]],Corrected=res[["dc"]])),col.regions=gcol(100),cuts=99,maxpixels=1e6,ylim=c(10,13),xlim=c(-7,-1)) |
|
| 60 |
} |
|
| 61 |
|
|
| 62 |
|
|
| 63 |
|
|
| 64 |
|
|
| 65 |
#### stuff below here is old junk..... |
|
| 66 |
|
|
| 67 |
## convert to netcdf, subset to mean/sd bands |
|
| 68 |
trans_ops=paste(" -co COMPRESS=DEFLATE -a_nodata 255 -stats -co FORMAT=NC4 -co ZLEVEL=9 -b 4 -b 2")
|
|
| 69 |
system(paste("gdal_translate -of netCDF ",trans_ops," ",ttif2," ",ncfile))
|
|
| 70 |
## file.remove(temptffile) |
|
| 71 |
system(paste("ncecat -O -u time ",ncfile," ",ncfile,sep=""))
|
|
| 72 |
## create temporary nc file with time information to append to CF data |
|
| 73 |
cat(paste("
|
|
| 74 |
netcdf time {
|
|
| 75 |
dimensions: |
|
| 76 |
time = 1 ; |
|
| 77 |
variables: |
|
| 78 |
int time(time) ; |
|
| 79 |
time:units = \"days since 2000-01-01 ",ifelse(s=="MOD09","10:30:00","13:30:00"),"\" ; |
|
| 80 |
time:calendar = \"gregorian\"; |
|
| 81 |
time:long_name = \"time of observation\"; |
|
| 82 |
data: |
|
| 83 |
time=",as.integer(date-as.Date("2000-01-01")),";
|
|
| 84 |
}"),file=paste(tempdir(),"/",date,"_time.cdl",sep="")) |
|
| 85 |
system(paste("ncgen -o ",tempdir(),"/",date,"_time.nc ",tempdir(),"/",date,"_time.cdl",sep=""))
|
|
| 86 |
## add time dimension to ncfile and compress (deflate) |
|
| 87 |
system(paste("ncks --fl_fmt=netcdf4 -L 9 -A ",tempdir(),"/",date,"_time.nc ",ncfile,sep=""))
|
|
| 88 |
## add other attributes |
|
| 89 |
system(paste("ncrename -v Band1,CF ",ncfile,sep=""))
|
|
| 90 |
system(paste("ncrename -v Band2,CFsd ",ncfile,sep=""))
|
|
| 91 |
## build correction factor explanation |
|
| 92 |
system(paste("ncatted ",
|
|
| 93 |
## CF Mean |
|
| 94 |
" -a units,CF,o,c,\"%\" ", |
|
| 95 |
" -a valid_range,CF,o,ub,\"0,100\" ", |
|
| 96 |
# " -a scale_factor,CF,o,b,\"0.1\" ", |
|
| 97 |
" -a _FillValue,CF,o,ub,\"255\" ", |
|
| 98 |
" -a missing_value,CF,o,ub,\"255\" ", |
|
| 99 |
" -a long_name,CF,o,c,\"Cloud Frequency (%)\" ", |
|
| 100 |
" -a correction_factor_description,CF,o,c,\"To account for variable observation frequency, CF in each pixel was adjusted by the proportion of days with at least one MODIS observation\" ", |
|
| 101 |
" -a correction_factor,CF,o,f,\"",round(modbeta1,4),"\" ", |
|
| 102 |
" -a NETCDF_VARNAME,CF,o,c,\"Cloud Frequency (%)\" ", |
|
| 103 |
## CF Standard Deviation |
|
| 104 |
" -a units,CFsd,o,c,\"SD\" ", |
|
| 105 |
" -a valid_range,CFsd,o,ub,\"0,200\" ", |
|
| 106 |
" -a scale_factor,CFsd,o,f,\"0.25\" ", |
|
| 107 |
" -a _FillValue,CFsd,o,ub,\"255\" ", |
|
| 108 |
" -a missing_value,CFsd,o,ub,\"255\" ", |
|
| 109 |
" -a long_name,CFsd,o,c,\"Cloud Frequency (%) Intra-month (2000-2013) Standard Deviation\" ", |
|
| 110 |
" -a NETCDF_VARNAME,CFsd,o,c,\"Cloud Frequency (%) Intra-month (2000-2013) Standard Deviation\" ", |
|
| 111 |
## global |
|
| 112 |
" -a title,global,o,c,\"Cloud Climatology from MOD09 Cloud Mask\" ", |
|
| 113 |
" -a institution,global,o,c,\"Jetz Lab, EEB, Yale University, New Haven, CT\" ", |
|
| 114 |
" -a source,global,o,c,\"Derived from MOD09GA Daily Data\" ", |
|
| 115 |
" -a comment,global,o,c,\"Developed by Adam M. Wilson (adam.wilson@yale.edu / http://adamwilson.us)\" ", |
|
| 116 |
ncfile,sep="")) |
|
| 117 |
|
|
| 118 |
## add the fillvalue attribute back (without changing the actual values) |
|
| 119 |
#system(paste("ncatted -a _FillValue,CF,o,b,-32768 ",ncfile,sep=""))
|
|
| 120 |
|
|
| 121 |
if(as.numeric(system(paste("cdo -s ntime ",ncfile),intern=T))<1) {
|
|
| 122 |
print(paste(ncfile," has no time, deleting")) |
|
| 123 |
file.remove(ncfile) |
|
| 124 |
} |
|
| 125 |
print(paste(basename(ncfile)," Finished")) |
|
| 126 |
|
|
| 127 |
|
|
| 128 |
} |
|
| 129 |
|
|
| 130 |
if(ramdisk) {
|
|
| 131 |
## unmount the ram disk |
|
| 132 |
system(paste("sudo umount ",tmpfs)
|
|
| 133 |
} |
|
| 134 |
|
|
| 135 |
|
|
| 136 |
### merge all the tiles to a single global composite |
|
| 137 |
#system(paste("ncdump -h ",list.files(tempdir,pattern="mod09.*.nc$",full=T)[10]))
|
|
| 138 |
file.remove("tmp/mod09_2000-01-15.nc")
|
|
| 139 |
system(paste("cdo -O mergetime -setrtomiss,-32768,-1 ",paste(list.files(tempdir,pattern="mod09.*.nc$",full=T),collapse=" ")," data/cloud_monthly.nc"))
|
|
| 140 |
|
|
| 141 |
# Overall mean |
|
| 142 |
system(paste("cdo -O timmean data/cloud_monthly.nc data/cloud_mean.nc"))
|
|
| 143 |
|
|
| 144 |
## Seasonal Means |
|
| 145 |
system(paste("cdo -f nc4c -O -yseasmean data/cloud_monthly.nc data/cloud_yseasmean.nc"))
|
|
| 146 |
system(paste("cdo -f nc4c -O -yseasstd data/cloud_monthly.nc data/cloud_yseasstd.nc"))
|
|
| 147 |
|
|
| 148 |
|
|
| 149 |
## standard deviation of mean monthly values give intra-annual variability |
|
| 150 |
system("cdo -f nc4c -z zip -chname,CF,CFsd -timstd data/cloud_ymonmean.nc data/cloud_std_intra.nc")
|
|
| 151 |
## mean of monthly standard deviations give inter-annual variability |
|
| 152 |
system("cdo -f nc4c -z zip -chname,CF,CFsd -timmean data/cloud_ymonstd.nc data/cloud_std_inter.nc")
|
|
| 153 |
|
|
| 154 |
|
|
| 155 |
|
|
| 156 |
|
|
| 157 |
## Daily animations |
|
| 158 |
regs=list( |
|
| 159 |
Venezuela=extent(c(-69,-59,0,7)), |
|
| 160 |
Cascades=extent(c(-122.8,-118,44.9,47)), |
|
| 161 |
Hawaii=extent(c(-156.5,-154,18.75,20.5)), |
|
| 162 |
Boliva=extent(c(-71,-63,-20,-15)), |
|
| 163 |
CFR=extent(c(17.75,22.5,-34.8,-32.6)), |
|
| 164 |
Madagascar=extent(c(46,52,-17,-12)) |
|
| 165 |
) |
|
| 166 |
|
|
| 167 |
r=1 |
|
| 168 |
|
|
| 169 |
system(paste("cdo -f nc4c -O inttime,2012-01-15,12:00:00,7day -sellonlatbox,",
|
|
| 170 |
paste(regs[[r]]@xmin,regs[[r]]@xmax,regs[[r]]@ymin,regs[[r]]@ymax,sep=","), |
|
| 171 |
" data/cloud_monthly.nc data/daily_",names(regs[r]),".nc",sep="")) |
|
| 172 |
|
|
| 173 |
|
|
| 174 |
|
|
| 175 |
|
|
| 176 |
### Long term summaries |
|
| 177 |
seasconc <- function(x,return.Pc=T,return.thetat=F) {
|
|
| 178 |
################################################################################################# |
|
| 179 |
## Precipitation Concentration function |
|
| 180 |
## This function calculates Precipitation Concentration based on Markham's (1970) technique as described in Schulze (1997) |
|
| 181 |
## South Africa Atlas of Agrohydology and Climatology - R E Schulze, M Maharaj, S D Lynch, B J Howe, and B Melvile-Thomson |
|
| 182 |
## Pages 37-38 |
|
| 183 |
################################################################################################# |
|
| 184 |
## x is a vector of precipitation quantities - the mean for each factor in "months" will be taken, |
|
| 185 |
## so it does not matter if the data are daily or monthly, as long as the "months" factor correctly |
|
| 186 |
## identifies them into 12 monthly bins, collapse indicates whether the data are already summarized as monthly means. |
|
| 187 |
################################################################################################# |
|
| 188 |
theta=seq(30,360,30)*(pi/180) # set up angles for each month & convert to radians |
|
| 189 |
if(sum(is.na(x))==12) { return(cbind(Pc=NA,thetat=NA)) ; stop}
|
|
| 190 |
if(return.Pc) {
|
|
| 191 |
rt=sqrt(sum(x * cos(theta))^2 + sum(x * sin(theta))^2) # the magnitude of the summation |
|
| 192 |
Pc=as.integer(round((rt/sum(x))*100))} |
|
| 193 |
if(return.thetat){
|
|
| 194 |
s1=sum(x*sin(theta),na.rm=T); s2=sum(x*cos(theta),na.rm=T) |
|
| 195 |
if(s1>=0 & s2>=0) {thetat=abs((180/pi)*(atan(sum(x*sin(theta),na.rm=T)/sum(x*cos(theta),na.rm=T))))}
|
|
| 196 |
if(s1>0 & s2<0) {thetat=180-abs((180/pi)*(atan(sum(x*sin(theta),na.rm=T)/sum(x*cos(theta),na.rm=T))))}
|
|
| 197 |
if(s1<0 & s2<0) {thetat=180+abs((180/pi)*(atan(sum(x*sin(theta),na.rm=T)/sum(x*cos(theta),na.rm=T))))}
|
|
| 198 |
if(s1<0 & s2>0) {thetat=360-abs((180/pi)*(atan(sum(x*sin(theta),na.rm=T)/sum(x*cos(theta),na.rm=T))))}
|
|
| 199 |
thetat=as.integer(round(thetat)) |
|
| 200 |
} |
|
| 201 |
if(return.thetat&return.Pc) return(c(conc=Pc,theta=thetat)) |
|
| 202 |
if(return.Pc) return(Pc) |
|
| 203 |
if(return.thetat) return(thetat) |
|
| 204 |
} |
|
| 205 |
|
|
| 206 |
|
|
| 207 |
|
|
| 208 |
## read in monthly dataset |
|
| 209 |
mod09=brick("data/cloud_ymonmean.nc",varname="CF")
|
|
| 210 |
plot(mod09[1]) |
|
| 211 |
|
|
| 212 |
mod09_seas=calc(mod09,seasconc,return.Pc=T,return.thetat=F,overwrite=T,filename="data/mod09_seas.nc",NAflag=255,datatype="INT1U") |
|
| 213 |
mod09_seas2=calc(mod09,seasconc,return.Pc=F,return.thetat=T,overwrite=T,filename="data/mod09_seas_theta.nc",datatype="INT1U") |
|
| 214 |
|
|
| 215 |
plot(mod09_seas) |
|
| 216 |
|
|
| 217 |
|
|
| climate/research/cloud/MOD09/biomes.R | ||
|---|---|---|
| 26 | 26 |
## add area and centroid to each polygon |
| 27 | 27 |
biome$areakm=do.call(c,mclapply(1:length(biome),function(i) {print(i); return(areaPolygon(biome[i,])/1000000)}))
|
| 28 | 28 |
biome@data[,c("lon","lat")]=coordinates(gCentroid(biome,byid=T))
|
| 29 |
## add numeric biome code for rasterization |
|
| 30 |
biome=biome[order(biome$realm,as.numeric(biome$biomeid)),] |
|
| 31 |
biome$icode=1:nrow(biome) |
|
| 32 |
## write it to disk as shapefile |
|
| 29 | 33 |
writeOGR(biome,"data/teow","biomes",driver="ESRI Shapefile",overwrite=T) |
| 30 | 34 |
} |
| 31 | 35 |
|
| 36 |
## rasterize biome shapefile to standard 1km grid |
|
| 37 |
ops=paste(" -ot Byte -at -a_nodata 255 -init 255 -a icode -te -180 -90 180 90 -tr 0.008333333333333 -0.008333333333333",
|
|
| 38 |
"-co BIGTIFF=yes -co COMPRESS=LZW -co PREDICTOR=1") |
|
| 39 |
system(paste("gdal_rasterize ",ops," -l biomes data/teow data/teow.tif"))
|
|
| 40 |
|
|
| 32 | 41 |
|
| 33 | 42 |
biome=readOGR("data/teow/","biomes")
|
| 34 | 43 |
projection(biome)="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs" |
| climate/research/cloud/MOD09/ee/ee_MCD09CF.js | ||
|---|---|---|
| 11 | 11 |
|
| 12 | 12 |
// Specify destination and run name |
| 13 | 13 |
var driveFolder="ee_mcd09cf"; |
| 14 |
var run="g2"
|
|
| 14 |
var run="g3"
|
|
| 15 | 15 |
|
| 16 | 16 |
// limit overall date range (only dates in this range will be included) |
| 17 | 17 |
var datestart=new Date("2000-01-01") // default time zone is UTC
|
| 18 |
var datestop=new Date("2014-02-28")
|
|
| 19 |
|
|
| 18 |
var datestop=new Date("2014-03-31")
|
|
| 20 | 19 |
|
| 21 | 20 |
// specify which months (within the date range above) to process |
| 22 | 21 |
var monthstart=1 |
| 23 | 22 |
var monthstop=12 |
| 24 | 23 |
|
| 24 |
// Sensors to process |
|
| 25 |
var mcols = ['MOD09GA','MYD09GA']; |
|
| 26 |
|
|
| 25 | 27 |
// specify what happens |
| 26 | 28 |
var verbose=false // print info about collections along the way (slows things down) |
| 27 |
var exportDrive=true // add exports to task window |
|
| 28 |
var drawmap=false // add image to map |
|
| 29 |
var drawmap=true // add image to map |
|
| 30 |
var test1=false // report on all images requested via dates above, but only add image below to map |
|
| 31 |
var exportDrive=!test1 // add exports to task window |
|
| 29 | 32 |
|
| 33 |
// set testing sensor and month, these only apply if test1==t above |
|
| 34 |
var testsensor='MYD09GA' |
|
| 35 |
var testmonth=1 |
|
| 30 | 36 |
|
| 31 | 37 |
// define regions |
| 32 | 38 |
var globe = '[[-180, -89.95], [-180, 89.95], [180, 89.95], [180, -89.95]]'; |
| ... | ... | |
| 44 | 50 |
var currentdate = new Date(); |
| 45 | 51 |
var date= currentdate.getFullYear()+''+("0" + (currentdate.getMonth() + 1)).slice(-2)+''+currentdate.getDate();
|
| 46 | 52 |
print(date) |
| 47 |
// identify start and stop years |
|
| 48 |
var yearstart=datestart.getUTCFullYear(); |
|
| 49 |
var yearstop=datestop.getUTCFullYear(); |
|
| 50 | 53 |
|
| 51 |
// get array of years to process |
|
| 52 |
var years=Array.apply(0, Array(yearstop-yearstart+1)).map(function (x, y) { return yearstart +y ; });
|
|
| 53 |
var nYears=years.length; |
|
| 54 |
if(verbose){print('Processing '+years)}
|
|
| 54 |
|
|
| 55 | 55 |
|
| 56 | 56 |
///////////////////////////////////////////////////////////////////////////// |
| 57 | 57 |
/// Functions |
| 58 | 58 |
|
| 59 |
//function to combine terra and aqua and limit by date
|
|
| 59 |
//function to select terra or aqua and subset by date
|
|
| 60 | 60 |
function getMCD09(modcol,datestart,datestop){
|
| 61 | 61 |
var getdates=ee.Filter.date(datestart,datestop); |
| 62 | 62 |
return(ee.ImageCollection(modcol).filter(getdates)); |
| ... | ... | |
| 136 | 136 |
select([0],['pObs'])}; // rename to nObs |
| 137 | 137 |
|
| 138 | 138 |
|
| 139 |
//function to return year+1 |
|
| 140 |
var yearplus=function (x, y) { return yearstart +y ; }
|
|
| 141 |
|
|
| 139 | 142 |
////////////////////////////////////////////////////////////////////////////////////////////////////// |
| 140 | 143 |
// Clip high-latitude nonsense |
| 141 | 144 |
// For an unknown reason there are pixels at high latitudes with nobs>0 in winter (dark) months. |
| ... | ... | |
| 145 | 148 |
var monthbox=function(month) {
|
| 146 | 149 |
// limits for each month (12 numbers correspond to jan-december limits) |
| 147 | 150 |
var ymin=[-90,-90,-90,-90,-70,-70,-70,-77,-77,-90,-90,-90]; |
| 148 |
var ymax=[ 77, 90, 90, 90, 90, 90, 90, 90, 90, 90, 77, 69];
|
|
| 151 |
var ymax=[ 75, 90, 90, 90, 90, 90, 90, 90, 90, 90, 77, 69];
|
|
| 149 | 152 |
// draw the polygon, use month-1 because month is 1:12 and array is indexed 0:11 |
| 150 | 153 |
var ind=month-1; |
| 151 | 154 |
var box = ee.Geometry.Rectangle(-180,ymin[ind],180,ymax[ind]); |
| ... | ... | |
| 158 | 161 |
//////////////////////////////////////////////////// |
| 159 | 162 |
// Start processing the data |
| 160 | 163 |
|
| 161 |
//var mcol='MOD09GA' |
|
| 162 |
//var tmonth=1 |
|
| 163 | 164 |
|
| 164 | 165 |
|
| 165 | 166 |
// loop over collections and months |
| 166 |
var mcols = ['MOD09GA','MYD09GA']; |
|
| 167 | 167 |
for (var c=0;c<mcols.length;c++) {
|
| 168 | 168 |
var mcol=mcols[c]; |
| 169 | 169 |
for (var tmonth = monthstart; tmonth <= monthstop; tmonth ++ ) {
|
| 170 | 170 |
|
| 171 |
// if test1 is true, use only test settings |
|
| 172 |
if(test1){
|
|
| 173 |
var mcol=testsensor |
|
| 174 |
var tmonth=testmonth |
|
| 175 |
} |
|
| 176 |
|
|
| 171 | 177 |
if(verbose){
|
| 172 | 178 |
print('Starting processing for:'+ mcol+' month '+tmonth);
|
| 173 | 179 |
} |
| 174 | 180 |
|
| 181 |
// identify start and stop years |
|
| 182 |
var yearstart=datestart.getUTCFullYear(); |
|
| 183 |
var yearstop=datestop.getUTCFullYear(); |
|
| 184 |
|
|
| 185 |
// update startdates based on terra/aqua start dates |
|
| 186 |
// otherwise missing month-years will cause problems when compiling the mean |
|
| 187 |
// Terra (MOD) startdate: February 24, 2000 |
|
| 188 |
if(mcol=='MOD09GA'&tmonth==1&yearstart==2000) yearstart=2001 |
|
| 189 |
/// Aqua (MYD) startdate: July 4, 2002 |
|
| 190 |
if(mcol=='MYD09GA'&tmonth<7&yearstart<=2002) yearstart=2003 |
|
| 191 |
if(mcol=='MYD09GA'&tmonth>=7&yearstart<=2002) yearstart=2002 |
|
| 192 |
|
|
| 193 |
// get array of years to process |
|
| 194 |
var years=Array.apply(0, Array(yearstop-yearstart+1)).map(yearplus); |
|
| 195 |
var nYears=years.length; |
|
| 196 |
if(verbose){print('Processing '+years)}
|
|
| 197 |
|
|
| 198 |
|
|
| 175 | 199 |
|
| 176 | 200 |
// build a combined M*D09GA collection limited by start and stop dates |
| 177 | 201 |
var MCD09all=getMCD09(mcol,datestart,datestop); |
| ... | ... | |
| 226 | 250 |
}); |
| 227 | 251 |
} |
| 228 | 252 |
|
| 253 |
if(test1) break; // if running testing, dont complete the loop |
|
| 254 |
|
|
| 229 | 255 |
} // close loop through months |
| 230 | 256 |
} // close loop through collections |
| 231 | 257 |
|
| ... | ... | |
| 239 | 265 |
addToMap(MCD09.select([2]),{min:0,max:15000,palette:palette}, "nObs");
|
| 240 | 266 |
addToMap(MCD09.select([3]),{min:0,max:10000,palette:palette}, "pObs");
|
| 241 | 267 |
// addToMap(trend.select('long-trend'), {min:-10, max:10, palette:palette}, 'trend');
|
| 242 |
} |
|
| 268 |
} |
|
| climate/research/cloud/MOD09/ee/ee_compile.R | ||
|---|---|---|
| 19 | 19 |
|
| 20 | 20 |
### Download files from google drive |
| 21 | 21 |
download=T |
| 22 |
if(download) system(paste("google docs get 2014032* ",datadir,"/mcd09ee",sep=""))
|
|
| 22 |
if(download) system(paste("google docs get 2014* ",datadir,"/mcd09ee",sep=""))
|
|
| 23 | 23 |
|
| 24 | 24 |
|
| 25 | 25 |
## Get list of available files |
| 26 |
df=data.frame(path=list.files(paste(datadir,"mcd09ee",sep="/"),pattern="*.tif$",full=T,recur=T),stringsAsFactors=F) |
|
| 26 |
version="g3" #which version of data from EE? |
|
| 27 |
df=data.frame(path=list.files(paste(datadir,"mcd09ee",sep="/"),pattern=paste(".*",version,".*.tif$",sep=""),full=T,recur=T),stringsAsFactors=F)
|
|
| 27 | 28 |
df[,c("month","sensor")]=do.call(rbind,strsplit(basename(df$path),"_|[.]|-"))[,c(5,4)]
|
| 28 | 29 |
df$date=as.Date(paste(2013,"_",df$month,"_15",sep=""),"%Y_%m_%d") |
| 29 | 30 |
|
| ... | ... | |
| 31 | 32 |
## use ramdisk? |
| 32 | 33 |
tmpfs="tmp/"#tempdir() |
| 33 | 34 |
|
| 34 |
|
|
| 35 | 35 |
ramdisk=F |
| 36 | 36 |
if(ramdisk) {
|
| 37 | 37 |
system("sudo mkdir -p /mnt/ram")
|
| ... | ... | |
| 45 | 45 |
|
| 46 | 46 |
rerun=T # set to true to recalculate all dates even if file already exists |
| 47 | 47 |
|
| 48 |
## define month-sensors to process |
|
| 49 |
jobs=unique(data.frame(month=df$month,sensor=df$sensor)) |
|
| 48 | 50 |
|
| 49 |
jobs=expand.grid(month=1:12,sensor=c("MOD09GA","MYD09GA"))
|
|
| 50 | 51 |
i=1 |
| 51 |
|
|
| 52 | 52 |
#jobs=jobs[jobs$sensor=="MYD09",] |
| 53 | 53 |
|
| 54 | 54 |
|
| ... | ... | |
| 96 | 96 |
|
| 97 | 97 |
|
| 98 | 98 |
|
| 99 |
|
|
| 100 |
|
|
| 101 |
## Create combined (MOD+MYD) uncorrected mean CF |
|
| 102 |
# foreach(i=1:12) %dopar% {
|
|
| 103 |
# ## get files |
|
| 104 |
# f=list.files("/mnt/data2/projects/cloud/mcd09tif",pattern=paste(".*[O|Y].*_",i,"[.]tif$",sep=""),full=T)
|
|
| 105 |
# ## Define output and check if it already exists |
|
| 106 |
# tmcd=paste("/mnt/data2/projects/cloud/mcd09tif/MCD09_",sprintf("%02d", i),".tif",sep="")
|
|
| 107 |
# ## check if output already exists |
|
| 108 |
# ops=paste("-t_srs 'EPSG:4326' -multi -srcnodata 255 -dstnodata 255 -r bilinear -te -180 -90 180 90 -tr 0.008333333333333 -0.008333333333333",
|
|
| 109 |
# "-co BIGTIFF=YES -ot Byte --config GDAL_CACHEMAX 500 -wm 500 -wo NUM_THREADS:10 -wo SOURCE_EXTRA=5") |
|
| 110 |
# system(paste("gdalwarp -overwrite -r average -co COMPRESS=LZW -co ZLEVEL=9 ",ops," ",paste(f,collapse=" ")," ",tmcd))
|
|
| 111 |
# ## update metadata |
|
| 112 |
# tags=c(paste("TIFFTAG_IMAGEDESCRIPTION='Monthly Cloud Frequency for 2000-2013 extracted from C5 MODIS MOD09GA and MYD09GA PGE11 internal cloud mask algorithm (embedded in state_1km bit 10).",
|
|
| 113 |
# "The daily cloud mask time series were summarized to mean cloud frequency (CF) by calculating the proportion of cloudy days. ", |
|
| 114 |
# "Band Descriptions: 1) Mean Monthly Cloud Frequency 2) Four Times the Standard Deviation of Mean Monthly Cloud Frequency", |
|
| 115 |
# " 3) Mean number of daily observations for each pixel 4) Proportion of days with at least one observation '"), |
|
| 116 |
# "TIFFTAG_DOCUMENTNAME='Collection 5 MCD09 Cloud Frequency'", |
|
| 117 |
# paste("TIFFTAG_DATETIME='2013",sprintf("%02d", i),"15'",sep=""),
|
|
| 118 |
# "TIFFTAG_ARTIST='Adam M. Wilson (adam.wilson@yale.edu)'") |
|
| 119 |
# system(paste("/usr/local/src/gdal-1.10.0/swig/python/scripts/gdal_edit.py ",tmcd," ",paste("-mo ",tags,sep="",collapse=" "),sep=""))
|
|
| 120 |
# writeLines(paste("Finished month",i))
|
|
| 121 |
# } |
|
| 122 |
|
|
| 123 |
|
|
| 124 |
### Perform bias correction |
|
| 125 |
foreach(i=1:nrow(jobs)) %dopar% {
|
|
| 126 |
## get month |
|
| 127 |
m=jobs$month[i] |
|
| 128 |
date=df$date[df$month==m][1] |
|
| 129 |
print(date) |
|
| 130 |
|
|
| 131 |
ttif1=paste(datadir,"/mcd09tif/",s,"_",sprintf("%02d", m),".tif",sep="")
|
|
| 132 |
ttif2=paste(tmpfs,"/",s,"_",m,"_wgs84.tif",sep="") |
|
| 133 |
ncfile=paste("data/mcd09nc/",s,"_",m,".nc",sep="")
|
|
| 134 |
|
|
| 135 |
## |
|
| 136 |
mod=stack(ttif1) |
|
| 137 |
names(mod)=c("cf","cfsd","nobs","pobs")
|
|
| 138 |
|
|
| 139 |
## set up processing chunks |
|
| 140 |
nrw=nrow(mod) |
|
| 141 |
nby=20 |
|
| 142 |
nrwg=seq(1,nrw,by=nby) |
|
| 143 |
writeLines(paste("Processing ",length(nrwg)," groups and",nrw,"lines"))
|
|
| 144 |
|
|
| 145 |
output=mclapply(nrwg,function(ti){
|
|
| 146 |
## Extract focal areas |
|
| 147 |
nr=51 |
|
| 148 |
nc=101 |
|
| 149 |
ngb=c(nr,nc) |
|
| 150 |
vals_cf=getValuesFocal(mod[[c("cf","pobs")]],ngb=ngb,row=ti,nrows=nby)
|
|
| 151 |
vals_obs=getValuesFocal(obs,ngb=ngb,row=ti,nrows=nby) |
|
| 152 |
## extract bias |
|
| 153 |
bias=raster(matrix(do.call(rbind,lapply(1:nrow(vals_cf),function(i) {
|
|
| 154 |
# if(i%in%round(seq(1,nrow(vals_cf),len=100))) print(i) |
|
| 155 |
tobs=vals_obs[i,] #vector of indices |
|
| 156 |
tval=vals_cf[i,] # vector of values |
|
| 157 |
lm1=lm(tval~tobs,na.rm=T) |
|
| 158 |
dif=round(predict(lm1)-predict(lm1,newdata=data.frame(tobs=median(tobs,na.rm=T)))) |
|
| 159 |
return(dif) |
|
| 160 |
})),nrow=nby,ncol=ncol(d),byrow=T)) # turn it back into a raster |
|
| 161 |
## update raster and write it |
|
| 162 |
extent(bias)=extent(d[ti:(ti+nby-1),1:ncol(d),drop=F]) |
|
| 163 |
projection(bias)=projection(d) |
|
| 164 |
NAvalue(bias) <- 255 |
|
| 165 |
writeRaster(bias,file=paste("data/bias/tiles/bias_",ti,".tif",sep=""),
|
|
| 166 |
format="GTiff",dataType="INT1U",overwrite=T,NAflag=255) #,options=c("COMPRESS=LZW","ZLEVEL=9")
|
|
| 167 |
print(ti) |
|
| 168 |
} |
|
| 169 |
) |
|
| 170 |
|
|
| 171 |
|
|
| 172 |
|
|
| 173 |
modpts=sampleRandom(cmod, size=10000, na.rm=TRUE, xy=T, sp=T) |
|
| 174 |
rm(cmod) #remove temporary raster to save space |
|
| 175 |
modpts=modpts[modpts$nobs>0,] #drop data from missing tiles |
|
| 176 |
### fit popbs correctionmodel (accounting for spatial variation in cf) |
|
| 177 |
modlm1=bam(cf~s(x,y)+pobs,data=modpts@data) |
|
| 178 |
summary(modlm1) |
|
| 179 |
modbeta1=coef(modlm1)["pobs"] |
|
| 180 |
|
|
| 181 |
writeLines(paste(date," slope:",round(modbeta1,4))) |
|
| 182 |
|
|
| 183 |
## mask no data regions (with less than 1 observation per day within that month) |
|
| 184 |
## use model above to correct for orbital artifacts |
|
| 185 |
biasf=function(cf,cfsd,nobs,pobs) {
|
|
| 186 |
## drop data in areass with nobs<1 |
|
| 187 |
drop=nobs<0|pobs<=50 |
|
| 188 |
cf[drop]=NA |
|
| 189 |
cfsd[drop]=NA |
|
| 190 |
nobs[drop]=NA |
|
| 191 |
pobs[drop]=NA |
|
| 192 |
bias=round((100-pobs)*modbeta1) |
|
| 193 |
cfc=cf+bias |
|
| 194 |
return(c(cf=cf,cfsd=cfsd,bias=bias,cfc=cfc,nobs=nobs,pobs=pobs))} |
|
| 195 |
|
|
| 196 |
# treg=extent(c(-1434564.00523,1784892.95369, 564861.173869, 1880991.3772)) |
|
| 197 |
# treg=extent(c(-2187230.72881, 4017838.07688, -339907.592509, 3589340.63805)) #all sahara |
|
| 198 |
|
|
| 199 |
mod2=overlay(cmod,fun=biasf,unstack=TRUE,filename=ttif1,format="GTiff", |
|
| 200 |
dataType="INT1U",overwrite=T,NAflag=255, options=c("COMPRESS=LZW", "BIGTIFF=YES"))
|
|
| 201 |
|
|
| 202 |
## warp to wgs84 |
|
| 203 |
ops=paste("-t_srs 'EPSG:4326' -multi -srcnodata 255 -dstnodata 255 -r bilinear -te -180 -90 180 90 -tr 0.008333333333333 -0.008333333333333",
|
|
| 204 |
"-co BIGTIFF=YES -ot Byte --config GDAL_CACHEMAX 500 -wm 500 -wo NUM_THREADS:10 -wo SOURCE_EXTRA=5") |
|
| 205 |
# ops=paste("-t_srs 'EPSG:4326' -srcnodata 255 -dstnodata 255 -multi -r bilinear -tr 0.008333333333333 -0.008333333333333",
|
|
| 206 |
# "-co BIGTIFF=YES -ot Byte --config GDAL_CACHEMAX 300000 -wm 300000 -wo NUM_THREADS:10 -wo SOURCE_EXTRA=5") |
|
| 207 |
system(paste("gdalwarp -overwrite ",ops," ",ttif1," ",ttif2))
|
|
| 99 |
#### stuff below here is old junk..... |
|
| 208 | 100 |
|
| 209 | 101 |
## convert to netcdf, subset to mean/sd bands |
| 210 | 102 |
trans_ops=paste(" -co COMPRESS=DEFLATE -a_nodata 255 -stats -co FORMAT=NC4 -co ZLEVEL=9 -b 4 -b 2")
|
Also available in: Unified diff
rearranging files to more linear progression through processing