Revision b74bf08c
Added by Adam Wilson almost 11 years ago
climate/research/cloud/MOD09/MOD09_Visualize.R | ||
---|---|---|
6 | 6 |
|
7 | 7 |
## read in global coasts for nice plotting |
8 | 8 |
library(maptools) |
9 |
library(doMC) |
|
9 | 10 |
|
11 |
registerDoMC(6) |
|
12 |
#### Evaluate MOD35 Cloud data |
|
13 |
#mc=brick("~/acrobates/adamw/projects/cloud/data/mod09.nc",varname="CF") |
|
14 |
f=list.files("/mnt/data2/projects/cloud/mcd09tif",pattern=paste(".*[O].*[.]tif$",sep=""),full=T) |
|
15 |
f=f[c(4:12,1:3)] |
|
16 |
mc=stack(f,bands=1) |
|
17 |
names(mc)=month.name |
|
18 |
|
|
19 |
#NAvalue(mc)=-1 |
|
10 | 20 |
#coast=spTransform(coast,CRS(projection(mod35))) |
11 | 21 |
land=readShapePoly("/mnt/data/jetzlab/Data/environ/global/gshhg/GSHHS_shp/c/GSHHS_c_L1.shp",force_ring=TRUE) |
12 | 22 |
projection(land)="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs" |
13 |
CP <- as(extent(-180, 180, -60, 84), "SpatialPolygons")
|
|
23 |
CP <- as(extent(-180, 180, -60, 66), "SpatialPolygons")
|
|
14 | 24 |
proj4string(CP) <- CRS(proj4string(land)) |
15 |
coast=as(land[land$area>50,],"SpatialLines") |
|
16 |
|
|
17 |
|
|
18 |
#### Evaluate MOD35 Cloud data |
|
19 |
mc=brick("~/acrobates/adamw/projects/cloud/data/mod09.nc",varname="CF") |
|
20 |
NAvalue(mc)=-1 |
|
21 | 25 |
|
22 |
cols=colorRampPalette(c("#000000","#00FF00","#FF0000"))#"black","blue","red")) |
|
23 |
for(i in 1:156){ |
|
24 |
png(paste("output/mod09_fullanimation_",i,".png",sep=""),width=2000,height=1000) |
|
25 |
print(i) |
|
26 |
r=mm[[i]] |
|
27 |
print(levelplot(r,col.regions=cols(100),at=seq(0,100,len=100),margin=F,maxpixels=1e6,ylim=c(-60,70),main=paste(names(mod09)[i])))+ |
|
28 |
layer(sp.lines(coast)) |
|
29 |
dev.off() |
|
30 |
} |
|
26 |
## World outline |
|
27 |
world=readOGR("/mnt/data/jetzlab/Data/environ/global/MODIS/modis_sinusoidal/","worldborder_sinusoidal") |
|
31 | 28 |
|
32 |
#### Evaluate MOD35 Cloud data
|
|
33 |
mmc=brick("~/acrobates/adamw/projects/cloud/data/cloud_ymonmean.nc",varname="CF")
|
|
34 |
names(mmc)=month.name
|
|
29 |
coast=crop(as(land[land$area>50,],"SpatialLines"),CP)
|
|
30 |
coast2=spTransform(coast,CRS(projection(mc)))
|
|
31 |
CP2=spTransform(CP,CRS(projection(mc)))
|
|
35 | 32 |
|
36 | 33 |
cols=colorRampPalette(c("#000000","#00FF00","#FF0000"))#"black","blue","red")) |
37 |
png("output/CF_Animation_%03d.png",width=5000,height=4000,res=600,pointsize=96,bg="white") |
|
38 |
for(i in 1:12){ |
|
34 |
foreach(i=1:12)%dopar% { |
|
39 | 35 |
print(i) |
40 |
r=mmc[[i]] |
|
41 |
print(levelplot(r,col.regions=cols(100),at=seq(1,100,len=100),margin=F,maxpixels=1e6,ylim=c(-60,73), |
|
36 |
r=mc[[i]] |
|
37 |
d1=(extent(coast2)@ymax-extent(coast2)@ymin)/(extent(coast2)@xmax-extent(coast2)@xmin) # get y:x ratio to set output |
|
38 |
png(paste("output/CF_mean_",sprintf("%02d", i),".png",sep=""),width=1920,height=round(1920*d1),res=300,pointsize=46,bg="white") |
|
39 |
print(levelplot(r,col.regions=cols(100),at=seq(1,100,len=100),margin=F,maxpixels=2.1e6,ylim=c(extent(coast2)@ymin,extent(coast2)@ymax), |
|
42 | 40 |
main=paste(month.name[i]),cex.main=3,scales=list(draw=F),cuts=99,ylab="",xlab="")+ |
43 |
layer(panel.polygon(x=c(-180,-180,180,180),y=c(-90,90,90,-90),col="black"),under=T)+ |
|
44 |
layer(sp.lines(coast,col="black"),under=F)) |
|
41 |
layer(panel.polygon(c(extent(world)@xmin,extent(world)@xmin,extent(world)@xmax,extent(world)@xmax),y=c(extent(world)@ymin,extent(world)@ymax,extent(world)@ymax,extent(world)@ymin),col="black"),under=T)+ |
|
42 |
layer(sp.lines(coast2,col="black"),under=F)+ |
|
43 |
layer(sp.lines(world,col="black"),under=F)) |
|
44 |
|
|
45 |
dev.off() |
|
45 | 46 |
} |
46 |
dev.off() |
|
47 | 47 |
|
48 | 48 |
|
49 |
|
|
50 |
|
|
51 |
#################################################################### |
|
52 |
### Regional Comparisons |
|
53 |
## Compare with worldclim and NPP |
|
54 |
wc=stack(as.list(paste("/mnt/data/jetzlab/Data/environ/global/worldclim/prec_",1:12,".bil",sep=""))) |
|
55 |
#wc_map=stack(as.list(paste("/mnt/data/jetzlab/Data/environ/global/worldclim/bio_12.bil",sep=""))) |
|
56 |
|
|
57 |
regs=list( |
|
58 |
Cascades=extent(c(-122.8,-118,44.9,47)), |
|
59 |
Hawaii=extent(c(-156.5,-154,18.75,20.5)), |
|
60 |
Boliva=extent(c(-71,-63,-20,-15)), |
|
61 |
Venezuela=extent(c(-69,-59,0,7)), |
|
62 |
CFR=extent(c(17.75,22.5,-34.8,-32.6)), |
|
63 |
Madagascar=extent(c(46,52,-17,-12)) |
|
64 |
#reg2=extent(c(-81,-70,-4,10)) |
|
65 |
) |
|
66 |
# convert to sinusoidal |
|
67 |
regs2=lapply(regs,function(r){ |
|
68 |
r2 = as(r, "SpatialPolygons") |
|
69 |
proj4string(r2) <- CRS(proj4string(land)) |
|
70 |
r3=spTransform(r2,CRS(projection(mc))) |
|
71 |
return(r3) |
|
72 |
}) |
|
73 |
|
|
74 |
|
|
75 |
## read in GEWEX 1-degree data |
|
76 |
gewex=brick("data/gewex/CA_PATMOSX_NOAA.nc",varname="a_CA") |
|
77 |
names(gewex)=month.name#"1-degree Cloud Frequency (PATMOS-x GEWEX AVHRR)" |
|
78 |
|
|
79 |
|
|
80 |
|
|
81 |
|
|
82 |
|
|
83 |
|
|
84 |
r="Venezuela" |
|
85 |
|
|
86 |
## print global map with box for region |
|
87 |
print(levelplot(mc[[1]],col.regions=cols(100),at=seq(1,100,len=100),margin=F,maxpixels=2.1e6,ylim=c(extent(coast2)@ymin,extent(coast2)@ymax), |
|
88 |
main=paste(month.name[i]),cex.main=3,scales=list(draw=F),cuts=99,ylab="",xlab="")+ |
|
89 |
layer(panel.polygon(c(extent(world)@xmin,extent(world)@xmin,extent(world)@xmax,extent(world)@xmax),y=c(extent(world)@ymin,extent(world)@ymax,extent(world)@ymax,extent(world)@ymin),col="black"),under=T)+ |
|
90 |
layer(sp.lines(coast2,col="black"),under=F)+ |
|
91 |
layer(sp.lines(world,col="black"),under=F)+ |
|
92 |
layer(sp.lines(regs[[r]],col="blue",lwd=2),under=F)) |
|
93 |
|
|
94 |
|
|
95 |
# ylab.right = "Cloud Frequency (%)",par.settings = list(layout.widths = list(axis.key.padding = 0.1,axis.left=0.6,ylab.right = 3,right.padding=2)), |
|
96 |
|
|
97 |
## crop the data |
|
98 |
tgewex=crop(gewex,regs[[r]]) |
|
99 |
tmap=crop(wc,regs[[r]]) |
|
100 |
tmc=projectRaster(crop(mc,regs2[[r]]),tmap) |
|
101 |
|
|
102 |
## set plotting parameters |
|
103 |
|
|
104 |
range(c(cellStats(tmc,range)),cellStats(tgewex,range)*100) |
|
105 |
range(cellStats(tmap,range)) |
|
106 |
|
|
107 |
pcols=colorRampPalette(c("#CCFFFF","#000066","#FF3300"))#"black","blue","red")) |
|
108 |
|
|
109 |
|
|
110 |
foreach(i=1:12)%dopar% { |
|
111 |
print(i) |
|
112 |
png(paste("output/CF_mean_",r,"_",sprintf("%02d", i),".png",sep=""),width=1920,height=1080,res=300,pointsize=46,bg="white") |
|
113 |
pars=list(layout.heights=list(key.bottom=2,key.top=1),layout.widths = list(axis.key.padding = 4,axis.left=1)) |
|
114 |
my.theme = trellis.par.get() |
|
115 |
my.theme$strip.background=list(col="transparent") |
|
116 |
trellis.par.set(my.theme) |
|
117 |
p1=levelplot(tmc[[i]],col.regions=cols(100),at=seq(20,100,len=99), |
|
118 |
colorkey=list(space="bottom",width=1,height=1.2,labels=list(labels=c(25,50,75,100),at=c(25,50,75,100)),just="left"), |
|
119 |
cuts=99,margin=F,max.pixels=1e6,par.settings = pars,main=paste(month.name[i]),cex.main=3,xlab="",ylab="",scales=list(x=list(at=c(-68,-64,-60)))) |
|
120 |
p2=levelplot(tgewex[[i]],col.regions=cols(100),at=seq(.20,1,len=99),cuts=99,margin=F,max.pixels=1e6, |
|
121 |
colorkey=F,#list(space="bottom",width=.01,height=0.00000001,labels=list(labels="",at="",cex=0.0001)), |
|
122 |
par.settings = pars) |
|
123 |
p3=levelplot(tmap[[i]],col.regions=pcols(100),cuts=100,at=seq(0,650,len=100),margin=F,maxpixels=1e6, |
|
124 |
colorkey=list(space="bottom",height=.5,width=1,labels=list(labels=c(0,300,600),at=c(0,300,600))),xlab="",ylab="",main=names(regs)[r],useRaster=T, |
|
125 |
par.settings = pars) |
|
126 |
p_all=c("MODIS Cloud (%)"=p1,"PATMOS-x Cloud (%)"=p2,"WorldClim Precip (mm)"=p3,x.same=T,y.same=T,merge.legends=T,layout=c(3,1)) |
|
127 |
print(p_all) |
|
128 |
dev.off() |
|
129 |
} |
climate/research/cloud/MOD09/ee/ee_compile.R | ||
---|---|---|
5 | 5 |
library(multicore) |
6 | 6 |
library(foreach) |
7 | 7 |
library(mgcv) |
8 |
registerDoMC(2) |
|
8 |
registerDoMC(12)
|
|
9 | 9 |
|
10 | 10 |
|
11 | 11 |
## start raster cluster |
12 |
#beginCluster(10)
|
|
12 |
#beginCluster(5)
|
|
13 | 13 |
|
14 | 14 |
setwd("~/acrobates/adamw/projects/cloud") |
15 | 15 |
|
... | ... | |
38 | 38 |
|
39 | 39 |
|
40 | 40 |
jobs=expand.grid(month=1:12,sensor=c("MOD09","MYD09")) |
41 |
## Loop over existing months to build composite netcdf files |
|
41 |
i=1 |
|
42 |
|
|
43 |
#jobs=jobs[jobs$sensor=="MYD09",] |
|
44 |
|
|
45 |
|
|
46 |
## Loop over data to mosaic tifs and remove problematic data at poles |
|
42 | 47 |
foreach(i=1:nrow(jobs)) %dopar% { |
43 | 48 |
## get month |
44 | 49 |
m=jobs$month[i] |
... | ... | |
47 | 52 |
## get sensor |
48 | 53 |
s=jobs$sensor[i] |
49 | 54 |
## Define output and check if it already exists |
50 |
tvrt=paste(tmpfs,"/",s,"_",m,".vrt",sep="") |
|
51 |
ttif1=paste(tmpfs,"/",s,"_",m,".tif",sep="") |
|
52 |
ttif2=paste(tmpfs,"/",s,"_",m,"_wgs84.tif",sep="") |
|
53 |
ncfile=paste("data/mcd09nc/",s,"_",m,".nc",sep="") |
|
55 |
tvrt=paste(tmpfs,"/",s,"_",sprintf("%02d", m),".vrt",sep="") |
|
56 |
ttif1=paste("data/mcd09tif/",s,"_",sprintf("%02d", m),".tif",sep="") |
|
54 | 57 |
|
55 | 58 |
## check if output already exists |
56 |
if(!rerun&file.exists(ncfile)) return(NA)
|
|
59 |
if(!rerun&file.exists(ttif1)) return(NA)
|
|
57 | 60 |
## build VRT to merge tiles |
58 | 61 |
system(paste("gdalbuildvrt ",tvrt," ",paste(df$path[df$month==m&df$sensor==s],collapse=" "))) |
59 |
|
|
60 |
##
|
|
62 |
|
|
63 |
## |
|
61 | 64 |
mod=stack(tvrt) |
62 | 65 |
names(mod)=c("cf","cfsd","nobs","pobs") |
63 |
################################# |
|
64 |
## correct for orbital artifacts |
|
65 | 66 |
|
66 |
## sample points to fit correction model |
|
67 |
reg=extent(-20016036,20016036, -4295789.46149, 4134293.61707) #banding region |
|
68 |
cmod=crop(mod,reg) |
|
69 |
names(cmod)=c("cf","cfsd","nobs","pobs") |
|
67 |
## mask no data regions (with less than 1 observation per day within that month) |
|
68 |
biasf=function(cf,cfsd,nobs,pobs) { |
|
69 |
## drop data in areass with nobs<1 or pobs<=50 or where sd=0 and cf=0 or 50 (some polar regions) |
|
70 |
drop=nobs<=0|pobs<=50|(cf==0&cfsd==0)|(cf==50&cfsd==0) |
|
71 |
cf[drop]=NA |
|
72 |
cfsd[drop]=NA |
|
73 |
return(c(cf=cf,cfsd=cfsd,nobs=nobs,pobs=pobs))} |
|
74 |
|
|
75 |
|
|
76 |
mod2=overlay(mod,fun=biasf,unstack=TRUE,filename=ttif1,format="GTiff", |
|
77 |
dataType="INT1U",overwrite=T,NAflag=255, options=c("COMPRESS=LZW","BIGTIFF=YES")) |
|
78 |
|
|
79 |
tags=c(paste("TIFFTAG_IMAGEDESCRIPTION='Monthly Cloud Frequency for 2000-2013 extracted from C5 MODIS ",s,"GA PGE11 internal cloud mask algorithm (embedded in state_1km bit 10).", |
|
80 |
"The daily cloud mask time series were summarized to mean cloud frequency (CF) by calculating the proportion of cloudy days. ", |
|
81 |
"Band Descriptions: 1) Mean Monthly Cloud Frequency 2) Four Times the Standard Deviation of Mean Monthly Cloud Frequency", |
|
82 |
" 3) Mean number of daily observations for each pixel 4) Proportion of days with at least one observation '"), |
|
83 |
"TIFFTAG_DOCUMENTNAME='Collection 5 ",s," Cloud Frequency'", |
|
84 |
paste("TIFFTAG_DATETIME='2013",sprintf("%02d", m),"15'",sep=""), |
|
85 |
"TIFFTAG_ARTIST='Adam M. Wilson (adam.wilson@yale.edu)'") |
|
86 |
system(paste("/usr/local/src/gdal-1.10.0/swig/python/scripts/gdal_edit.py ",ttif1," ",paste("-mo ",tags,sep="",collapse=" "),sep="")) |
|
87 |
|
|
88 |
rm(mod,mod2) |
|
89 |
writeLines(paste("Month:",m," Sensor:",s," Finished")) |
|
90 |
|
|
91 |
} |
|
92 |
|
|
93 |
|
|
94 |
## Create combined (MOD+MYD) uncorrected mean CF |
|
95 |
# foreach(i=1:12) %dopar% { |
|
96 |
# ## get files |
|
97 |
# f=list.files("/mnt/data2/projects/cloud/mcd09tif",pattern=paste(".*[O|Y].*_",i,"[.]tif$",sep=""),full=T) |
|
98 |
# ## Define output and check if it already exists |
|
99 |
# tmcd=paste("/mnt/data2/projects/cloud/mcd09tif/MCD09_",sprintf("%02d", i),".tif",sep="") |
|
100 |
# ## check if output already exists |
|
101 |
# ops=paste("-t_srs 'EPSG:4326' -multi -srcnodata 255 -dstnodata 255 -r bilinear -te -180 -90 180 90 -tr 0.008333333333333 -0.008333333333333", |
|
102 |
# "-co BIGTIFF=YES -ot Byte --config GDAL_CACHEMAX 500 -wm 500 -wo NUM_THREADS:10 -wo SOURCE_EXTRA=5") |
|
103 |
# system(paste("gdalwarp -overwrite -r average -co COMPRESS=LZW -co ZLEVEL=9 ",ops," ",paste(f,collapse=" ")," ",tmcd)) |
|
104 |
# ## update metadata |
|
105 |
# 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).", |
|
106 |
# "The daily cloud mask time series were summarized to mean cloud frequency (CF) by calculating the proportion of cloudy days. ", |
|
107 |
# "Band Descriptions: 1) Mean Monthly Cloud Frequency 2) Four Times the Standard Deviation of Mean Monthly Cloud Frequency", |
|
108 |
# " 3) Mean number of daily observations for each pixel 4) Proportion of days with at least one observation '"), |
|
109 |
# "TIFFTAG_DOCUMENTNAME='Collection 5 MCD09 Cloud Frequency'", |
|
110 |
# paste("TIFFTAG_DATETIME='2013",sprintf("%02d", i),"15'",sep=""), |
|
111 |
# "TIFFTAG_ARTIST='Adam M. Wilson (adam.wilson@yale.edu)'") |
|
112 |
# system(paste("/usr/local/src/gdal-1.10.0/swig/python/scripts/gdal_edit.py ",tmcd," ",paste("-mo ",tags,sep="",collapse=" "),sep="")) |
|
113 |
# writeLines(paste("Finished month",i)) |
|
114 |
# } |
|
115 |
|
|
116 |
|
|
117 |
### Perform bias correction |
|
118 |
foreach(i=1:nrow(jobs)) %dopar% { |
|
119 |
## get month |
|
120 |
m=jobs$month[i] |
|
121 |
date=df$date[df$month==m][1] |
|
122 |
print(date) |
|
123 |
|
|
124 |
ttif2=paste(tmpfs,"/",s,"_",m,"_wgs84.tif",sep="") |
|
125 |
ncfile=paste("data/mcd09nc/",s,"_",m,".nc",sep="") |
|
126 |
|
|
70 | 127 |
modpts=sampleRandom(cmod, size=10000, na.rm=TRUE, xy=T, sp=T) |
128 |
rm(cmod) #remove temporary raster to save space |
|
71 | 129 |
modpts=modpts[modpts$nobs>0,] #drop data from missing tiles |
72 | 130 |
### fit popbs correctionmodel (accounting for spatial variation in cf) |
73 | 131 |
modlm1=bam(cf~s(x,y)+pobs,data=modpts@data) |
74 | 132 |
summary(modlm1) |
75 | 133 |
modbeta1=coef(modlm1)["pobs"] |
76 | 134 |
|
77 |
|
|
78 |
#modlm2=bam(cf~s(x,y),data=modpts@data) |
|
79 |
#resid=residuals(modlm2)#,newdata=data.frame(x=modpts$x,y=modpts$y,pobs=100)) |
|
80 |
#pred=predict(modlm1,newdata=data.frame(x=modpts$x,y=modpts$y,pobs=100)) |
|
81 |
#plot(resid~pobs,data=modpts@data);abline(lm(resid~modpts$pobs)) |
|
82 |
#plot(modlm1);points(modpts) |
|
83 |
|
|
84 |
#writeLines(paste(date," slope:",round(modbeta1,4))) |
|
135 |
writeLines(paste(date," slope:",round(modbeta1,4))) |
|
136 |
## Smooth data to remove large-grain variability |
|
137 |
# system(paste("pkfilter -f mean -dx 99 -dy 99 -ot Byte -i ",df$path[df$month==m&df$sensor==s][1]," -o ",ttif1)) |
|
85 | 138 |
|
86 | 139 |
## mask no data regions (with less than 1 observation per day within that month) |
87 | 140 |
## use model above to correct for orbital artifacts |
... | ... | |
101 | 154 |
|
102 | 155 |
mod2=overlay(cmod,fun=biasf,unstack=TRUE,filename=ttif1,format="GTiff", |
103 | 156 |
dataType="INT1U",overwrite=T,NAflag=255, options=c("COMPRESS=LZW", "BIGTIFF=YES")) |
104 |
|
|
105 |
# system(paste("pksetmask -i ",modvrt," -m ",nObs, |
|
106 |
# " --operator='<' --msknodata 1 --nodata 255 --operator='>' --msknodata 10 --nodata 10 -o ",modtif1,sep="")) |
|
107 |
|
|
157 |
|
|
108 | 158 |
## warp to wgs84 |
109 | 159 |
ops=paste("-t_srs 'EPSG:4326' -multi -srcnodata 255 -dstnodata 255 -r bilinear -te -180 -90 180 90 -tr 0.008333333333333 -0.008333333333333", |
110 | 160 |
"-co BIGTIFF=YES -ot Byte --config GDAL_CACHEMAX 500 -wm 500 -wo NUM_THREADS:10 -wo SOURCE_EXTRA=5") |
climate/research/cloud/MOD09/ee_biastest.R | ||
---|---|---|
16 | 16 |
## set raster options |
17 | 17 |
rasterOptions(maxmemory=1e8) |
18 | 18 |
|
19 |
## get CF data |
|
20 |
## Get list of available files |
|
21 |
#df=data.frame(path=list.files("data/mod09cloud",pattern="*.tif$",full=T,recur=T),stringsAsFactors=F) |
|
22 |
#df[,c("month")]=do.call(rbind,strsplit(basename(df$path),"_|[.]|-"))[,c(6)] |
|
23 |
#df$date=as.Date(paste(2013,"_",df$month,"_15",sep=""),"%Y_%m_%d") |
|
24 |
#df=df[order(df$date),] |
|
19 |
#reg=extent(c(-1434564.00523,1784892.95369, 564861.173869, 1880991.3772)) |
|
20 |
#reg=extent(c(-2187230.72881, 4017838.07688, -339907.592509, 3589340.63805)) #all sahara |
|
21 |
#reg=extent(c(-12020769.9608, 10201058.231, -3682105.25271, 3649806.08382)) #large equitorial |
|
25 | 22 |
|
26 |
#d=stack(df$path,bands=1) |
|
27 |
#names(d)=df$month |
|
28 | 23 |
|
29 |
reg=extent(c(-1434564.00523,1784892.95369, 564861.173869, 1880991.3772)) |
|
30 |
reg=extent(c(-2187230.72881, 4017838.07688, -339907.592509, 3589340.63805)) #all sahara |
|
31 |
reg=extent(c(-12020769.9608, 10201058.231, -3682105.25271, 3649806.08382)) #large equitorial |
|
32 |
|
|
33 |
#reg=extent(c(-151270.082307, 557321.958761, 1246351.8356, 1733123.75947)) |
|
34 |
|
|
35 |
d=stack("data/mcd09/2014113_combined_20002013_MOD09_1-img-0000000000-0000000000.tif") |
|
36 |
|
|
37 |
|
|
38 |
names(d)=c("cf","cfsd","nobs","pobs") |
|
39 |
#cd=coordinates(d) |
|
40 |
#cd$x=as.integer(cd[,"x"]); cd$y=as.integer(cd[,"y"]) |
|
41 |
#cdr=rasterFromXYZ(cbind(cd,cd)) |
|
42 |
#names(cdr) |
|
43 |
#d=stack(d,cdr) |
|
44 |
#obs=crop(raster("data/mcd09/2014113_combined_20002013_MOD09_1-img-0000000000-0000000000.tif",band=4),reg) |
|
45 |
#levelplot(stack(obs,d)) |
|
46 |
|
|
47 |
#plot(obs,d) |
|
48 |
|
|
49 |
#pts=sampleStratified(d[["pobs"]], size=10000, exp=10, na.rm=TRUE, xy=FALSE, ext=NULL, sp=FALSE) |
|
50 |
pts=sampleRandom(d, size=100000, na.rm=TRUE, xy=T, ext=NULL, sp=T) |
|
51 |
pts=pts[pts$nobs>0,] #drop data from missing tiles |
|
52 |
|
|
53 |
#dt=data.frame(cf=as.numeric(values(d[["cf"]])),obs=as.numeric(values(d[["pobs"]]))) |
|
54 |
#dt[,c("x","y")]=coordinates(d) |
|
55 |
#dt=dt[dt$obs>0,] #drop data from missing tiles |
|
56 |
|
|
57 |
#s=1:nrow(dt) |
|
58 |
#s=sample(1:nrow(dt),10000) |
|
59 |
#s=strata(dt, stratanames=dt$obs, size=1000) |
|
60 |
|
|
61 |
lm1=bam(cf~s(x,y)+pobs+nobs,data=pts@data) |
|
62 |
summary(lm1) |
|
63 |
beta1=coef(lm1)["pobs"]; beta1 |
|
64 |
beta2=coef(lm1)["nobs"]; beta2 |
|
65 |
|
|
66 |
#d2=d |
|
67 |
#d2[["pobs"]]=100 |
|
68 |
|
|
69 |
#pred1=raster::predict(d,lm1,file=paste("data/bias/pred1.tif",sep=""),format="GTiff",dataType="INT1U",overwrite=T,NAflag=255)#lm1=bam(as.vector(d[["cf"]]) ~ as.vector(d[["pobs"]])) |
|
70 |
#pred2=raster::predict(d2,lm1,file=paste("data/bias/pred2.tif",sep=""),format="GTiff",dataType="INT1U",overwrite=T,NAflag=255)#lm1=bam(as.vector(d[["cf"]]) ~ as.vector(d[["pobs"]])) |
|
71 |
#pred3=overlay(pred1,pred2,function(x,y) x-y,file="data/bias/pred3.tif",overwrite=T) |
|
72 |
#biasc=overlay(d[["cf"]], pred3,fun=sum,file="data/bias/biasc.tif",overwrite=T) |
|
73 |
|
|
74 |
|
|
75 |
getbias=function(cf,cfsd,nobs,pobs) return((100-pobs)*beta1+(4-nobs)*beta2) |
|
76 |
getbias(cf=c(50,40,30),c(1,1,1),c(3,4,5),c(82,87,87)) |
|
77 |
bias=overlay(d,fun=getbias, unstack=TRUE,file=paste("data/bias/bias.tif",sep=""),format="GTiff",dataType="INT1U",overwrite=T,NAflag=255) |
|
78 |
dc=overlay(d[["cf"]],bias,fun=function(x,y) x+y, |
|
79 |
file=paste("data/bias/biasc.tif",sep=""), |
|
80 |
format="GTiff",dataType="INT1U",overwrite=T,NAflag=255) #,options=c("COMPRESS=LZW","ZLEVEL=9")) |
|
81 |
#levelplot(stack(d,dc)) |
|
82 |
|
|
83 |
|
|
84 |
library(multicore) |
|
85 |
|
|
86 |
Mode <- function(x) { |
|
87 |
ux <- unique(x) |
|
88 |
ux[which.max(tabulate(match(x, ux)))] |
|
89 |
} |
|
90 |
#rasterOptions(maxmemory)#=50GB) |
|
91 | 24 |
## set up processing chunks |
92 |
nrw=nrow(d) |
|
25 |
nrw=nrow(mod)
|
|
93 | 26 |
nby=20 |
94 | 27 |
nrwg=seq(1,nrw,by=nby) |
95 | 28 |
writeLines(paste("Processing ",length(nrwg)," groups and",nrw,"lines")) |
96 | 29 |
|
97 | 30 |
## Parallel loop to conduct moving window analysis and quantify pixels with significant shifts across pp or lulc boundaries |
98 | 31 |
output=mclapply(nrwg,function(ti){ |
99 |
## Extract focal areas |
|
100 |
ngb=c(51,101) |
|
101 |
vals_cf=getValuesFocal(d,ngb=ngb,row=ti,nrows=nby) |
|
32 |
## Extract focal areas |
|
33 |
nr=51 |
|
34 |
nc=101 |
|
35 |
ngb=c(nr,nc) |
|
36 |
vals_cf=getValuesFocal(mod[[c("cf","pobs")]],ngb=ngb,row=ti,nrows=nby) |
|
102 | 37 |
vals_obs=getValuesFocal(obs,ngb=ngb,row=ti,nrows=nby) |
103 | 38 |
## extract bias |
104 | 39 |
bias=raster(matrix(do.call(rbind,lapply(1:nrow(vals_cf),function(i) { |
Also available in: Unified diff
adding uncorrected MCD output in wgs84