Revision f0375bec
Added by Adam M. Wilson almost 11 years ago
| climate/research/cloud/MOD09/ee_biastest.R | ||
|---|---|---|
| 12 | 12 |
|
| 13 | 13 |
setwd("~/acrobates/adamw/projects/cloud")
|
| 14 | 14 |
|
| 15 |
|
|
| 15 | 16 |
## set raster options |
| 16 | 17 |
rasterOptions(maxmemory=1e8) |
| 17 | 18 |
|
| ... | ... | |
| 31 | 32 |
#reg=extent(c(-151270.082307, 557321.958761, 1246351.8356, 1733123.75947)) |
| 32 | 33 |
|
| 33 | 34 |
d=crop(stack("data/mcd09/2014113_combined_20002013_MOD09_1-img-0000000000-0000000000.tif"),reg)
|
| 35 |
d=crop(stack("/Users/adamw/GoogleDrive/EarthEngineOutput/ee_combined/2014113_combined_20002013_MOD09_1-img-0000000000-0000000000.tif"),reg)
|
|
| 36 |
|
|
| 34 | 37 |
names(d)=c("cf","cfsd","nobs","pobs")
|
| 38 |
cd=coordinates(d) |
|
| 39 |
cdr=rasterFromXYZ(cbind(cd,cd)) |
|
| 40 |
names(cdr) |
|
| 41 |
d=stack(d,cdr) |
|
| 35 | 42 |
#obs=crop(raster("data/mcd09/2014113_combined_20002013_MOD09_1-img-0000000000-0000000000.tif",band=4),reg)
|
| 36 | 43 |
#levelplot(stack(obs,d)) |
| 37 | 44 |
|
| ... | ... | |
| 51 | 58 |
|
| 52 | 59 |
lm1=bam(cf~s(x,y,pobs),data=pts@data) |
| 53 | 60 |
summary(lm1) |
| 54 |
beta=coef(lm1)["pobs"]; beta |
|
| 55 |
|
|
| 56 |
#lm1=bam(as.vector(d[["cf"]]) ~ as.vector(d[["pobs"]])) |
|
| 57 |
#fun <- function(x) { bam(x[["cf"]] ~ x[["pobs"]])$coefficients["obs"] }
|
|
| 58 |
#bias <- calc(d, fun) |
|
| 59 |
|
|
| 60 |
#pred=interpolate(d,lm1) |
|
| 61 |
|
|
| 62 |
getbias=function(x) return((100-x)*beta) |
|
| 63 |
|
|
| 64 |
|
|
| 65 |
getbias(beta,87) |
|
| 66 |
|
|
| 67 |
|
|
| 68 |
bias=calc(d[["pobs"]],getbias,file=paste("data/bias/bias.tif",sep=""),format="GTiff",dataType="INT1U",overwrite=T,NAflag=255)
|
|
| 69 |
|
|
| 70 |
dc=overlay(d[["cf"]],bias,fun=function(x,y) x+y, |
|
| 71 |
file=paste("data/bias/biasc.tif",sep=""),
|
|
| 72 |
format="GTiff",dataType="INT1U",overwrite=T,NAflag=255) #,options=c("COMPRESS=LZW","ZLEVEL=9"))
|
|
| 73 |
|
|
| 74 |
|
|
| 61 |
#beta=coef(lm1)["pobs"]; beta |
|
| 62 |
d2=d |
|
| 63 |
d2[["pobs"]]=100 |
|
| 64 |
|
|
| 65 |
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"]]))
|
|
| 66 |
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"]]))
|
|
| 67 |
pred3=overlay(pred1,pred2,function(x,y) x-y,file="data/bias/pred3.tif",overwrite=T) |
|
| 68 |
biasc=overlay(d[["cf"]], pred3,fun=sum,file="data/bias/biasc.tif",overwrite=T) |
|
| 69 |
|
|
| 70 |
|
|
| 71 |
#getbias=function(x) return((100-x)*beta) |
|
| 72 |
#getbias(beta,87) |
|
| 73 |
#bias=calc(d[["pobs"]],getbias,file=paste("data/bias/bias.tif",sep=""),format="GTiff",dataType="INT1U",overwrite=T,NAflag=255)
|
|
| 74 |
#dc=overlay(d[["cf"]],bias,fun=function(x,y) x+y, |
|
| 75 |
# file=paste("data/bias/biasc.tif",sep=""),
|
|
| 76 |
# format="GTiff",dataType="INT1U",overwrite=T,NAflag=255) #,options=c("COMPRESS=LZW","ZLEVEL=9"))
|
|
| 75 | 77 |
#levelplot(stack(d,dc)) |
| 76 | 78 |
|
| 77 | 79 |
|
Also available in: Unified diff
attempting gam to remove bias