Revision 598244e1
Added by Adam Wilson over 11 years ago
| climate/procedures/MOD35C5_Evaluation.r | ||
|---|---|---|
| 22 | 22 |
|
| 23 | 23 |
## mod35C6 annual |
| 24 | 24 |
if(!file.exists("data/MOD35C6_2009.tif")){
|
| 25 |
system("/usr/local/gdal-1.10.0/bin/gdalbuildvrt -a_srs '+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs' -sd 1 -b 1 data/MOD35C6.vrt /home/adamw/acrobates/adamw/projects/interp/data/modis/mod35/summary/*mean.nc ")
|
|
| 26 |
system("align.sh data/MOD35C6.vrt data/MOD09_2009.tif data/MOD35C6_2009.kmz")
|
|
| 25 |
system("/usr/local/gdal-1.10.0/bin/gdalbuildvrt -a_srs '+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs' -sd 1 -b 1 data/MOD35C6.vrt `find /home/adamw/acrobates/adamw/projects/interp/data/modis/mod35/summary/ -name '*h[1]*_mean.nc'` ")
|
|
| 26 |
system("align.sh data/MOD35C6.vrt data/MOD09_2009.tif data/MOD35C6_2009.tif")
|
|
| 27 | 27 |
} |
| 28 |
mod35c6=raster("data/MOD35C6_2009_v1.tif")
|
|
| 28 |
mod35c6=raster("data/MOD35C6_2009.tif")
|
|
| 29 | 29 |
names(mod35c6)="C6MOD35CF" |
| 30 | 30 |
NAvalue(mod35c6)=255 |
| 31 | 31 |
|
| 32 |
|
|
| 33 | 32 |
## landcover |
| 34 | 33 |
if(!file.exists("data/MCD12Q1_IGBP_2009_051_wgs84_1km.tif")){
|
| 35 | 34 |
system(paste("/usr/local/gdal-1.10.0/bin/gdalwarp -tr 0.008983153 0.008983153 -r mode -ot Byte -co \"COMPRESS=LZW\"",
|
| ... | ... | |
| 226 | 225 |
## exploring various ways to compare cloud products |
| 227 | 226 |
t1=trd1[[1]] |
| 228 | 227 |
dif_p=calc(trd1[[1]], function(x) (x[1]-x[3])/(1-x[1])) |
| 229 |
edge=edge(subset(t1,"MCD12Q1"),classes=T,type="inner") |
|
| 228 |
edge=calc(edge(subset(t1,"MCD12Q1"),classes=T,type="inner"),function(x) ifelse(x==1,1,NA)) |
|
| 229 |
edgeb=buffer(edge,width=5000) |
|
| 230 |
edgeb=calc(edgeb,function(x) ifelse(is.na(x),0,1)) |
|
| 231 |
|
|
| 230 | 232 |
names(edge)="edge" |
| 231 |
td1=as.data.frame(stack(t1,edge)) |
|
| 233 |
names(edgeb)="edgeb" |
|
| 234 |
|
|
| 235 |
td1=as.data.frame(stack(t1,edge,edgeb)) |
|
| 232 | 236 |
|
| 233 | 237 |
cor(td1$MOD17,td1$C5MOD35,use="complete",method="spearman") |
| 234 |
cor(td1$MOD17[td1$edge==1],td1$C5MOD35[td1$edge==1],use="complete",method="spearman") |
|
| 235 |
cor(td1,use="complete",method="spearman") |
|
| 238 |
cor(td1$MOD17[td1$edgeb==1],td1$C5MOD35[td1$edgeb==1],use="complete",method="spearman") |
|
| 239 |
round(cor(td1,use="complete",method="spearman"),2) |
|
| 240 |
## tests |
|
| 241 |
cor.test(td1$MOD17,td1$C5MOD35,use="complete",method="spearman",alternative="two.sided") |
|
| 242 |
cor.test(td1$MOD17,td1$C6MOD35,use="complete",method="spearman",alternative="two.sided") |
|
| 243 |
cor.test(td1$MOD17,td1$C5MOD35,use="complete",method="spearman",alternative="two.sided") |
|
| 244 |
|
|
| 245 |
cor.test( |
|
| 236 | 246 |
splom(t1) |
| 237 | 247 |
plot(mod17,mod17qc) |
| 238 |
xyplot(MOD17~C5MOD35CF|edge,data=td1) |
|
| 248 |
xyplot(MOD17~C5MOD35CF|edgeb,data=td1) |
|
| 249 |
bwplot(MCD12Q1~C5MOD35CF|edgeb,data=td1) |
|
| 250 |
|
|
| 239 | 251 |
plot(dif_p) |
| 240 | 252 |
|
| 241 | 253 |
#rondonia=trd[trd$trans=="Brazil",] |
| ... | ... | |
| 251 | 263 |
coast=map2SpatialLines(map("world", interior=FALSE, plot=FALSE),proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
|
| 252 | 264 |
|
| 253 | 265 |
g1=levelplot(stack(mod35c5,mod09),xlab=" ",scales=list(x=list(draw=F),y=list(alternating=1)),col.regions=cols,at=at)+layer(sp.polygons(bbs[1:4],lwd=2))+layer(sp.lines(coast,lwd=.5)) |
| 266 |
|
|
| 254 | 267 |
g2=levelplot(dif_c5_09,col.regions=bgrayr(100),at=seq(-70,70,len=100),margin=F,ylab=" ",colorkey=list("right"))+layer(sp.polygons(bbs[1:4],lwd=2))+layer(sp.lines(coast,lwd=.5))
|
| 255 |
trellis.par.set(background=list(fill="white"),panel.background=list(fill="white")) |
|
| 256 |
g3=histogram(dif_c5_09,bg="white",col="black",border=NA,scales=list(x=list(at=c(-50,0,50)),y=list(draw=F),cex=1))+layer(panel.abline(v=0,col="red",lwd=2)) |
|
| 268 |
g2$strip=strip.custom(var.name="Difference (C5MOD35-MOD09)",style=1,strip.names=T,strip.levels=F) #update strip text |
|
| 269 |
bg <- trellis.par.get("panel.background")
|
|
| 270 |
bg$col <- "white" |
|
| 271 |
trellis.par.set("panel.background",bg)
|
|
| 272 |
g3=histogram(dif_c5_09,col="black",border=NA,scales=list(x=list(at=c(-50,0,50)),y=list(draw=F),cex=1))+layer(panel.abline(v=0,col="red",lwd=2)) |
|
| 257 | 273 |
|
| 258 | 274 |
### regional plots |
| 259 | 275 |
p1=useOuterStrips(levelplot(value~x*y|variable+trans,data=trd[!trd$variable%in%c("MCD12Q1","MOD35pp"),],asp=1,scales=list(draw=F,rot=0,relation="free"),
|
| ... | ... | |
| 327 | 343 |
rectangles=list(border=NA,col=c(NA,IGBP$col[sort(unique(transd$value[transd$variable=="MCD12Q1"]+1))])), |
| 328 | 344 |
text=list(c("MCD12Q1 IGBP Land Cover",IGBP$class[sort(unique(transd$value[transd$variable=="MCD12Q1"]+1))])))))),
|
| 329 | 345 |
strip = strip.custom(par.strip.text=list(cex=.75))) |
| 330 |
print(p4) |
|
| 346 |
#print(p4)
|
|
| 331 | 347 |
|
| 332 | 348 |
#trdw=cast(trd,trans+x+y~variable,value="value") |
| 333 | 349 |
#transdw=cast(transd,transect+dist~variable,value="value") |
| ... | ... | |
| 342 | 358 |
#CairoPNG("output/mod35compare_%d.png",units="in", width=11,height=8.5,pointsize=4000,dpi=1200,antialias="subpixel")
|
| 343 | 359 |
### Global Comparison |
| 344 | 360 |
print(g1) |
| 345 |
print(g1,position=c(0,.33,1,1),more=T) |
|
| 346 |
print(g2,position=c(0,0,1,0.394),more=T) |
|
| 347 |
print(g3,position=c(0.31,0.06,.42,0.27),more=F) |
|
| 361 |
print(g1,position=c(0,.35,1,1),more=T) |
|
| 362 |
print(g2,position=c(0,0,1,0.415),more=F) |
|
| 363 |
#print(g3,position=c(0.31,0.06,.42,0.27),more=F) |
|
| 364 |
|
|
| 348 | 365 |
### MOD35 Desert Processing path |
| 349 | 366 |
levelplot(pp,asp=1,scales=list(draw=T,rot=0),maxpixels=1e6, |
| 350 | 367 |
at=c(-1:3),col.regions=c("blue","cyan","tan","darkgreen"),margin=F,
|
| ... | ... | |
| 358 | 375 |
dev.off() |
| 359 | 376 |
|
| 360 | 377 |
### summary stats for paper |
| 361 |
td=cast(transect+loc+dist~variable,value="value",data=transd)
|
|
| 362 |
td2=melt.data.frame(td,id.vars=c("transect","dist","loc","MOD35pp","MCD12Q1"))
|
|
| 378 |
td=cast(transect+loc+dist~variable,value="value",data=transd) |
|
| 379 |
td2=melt.data.frame(td,id.vars=c("transect","dist","loc","MOD35pp","MCD12Q1"))
|
|
| 363 | 380 |
|
| 364 | 381 |
## function to prettyprint mean/sd's |
| 365 | 382 |
msd= function(x) paste(round(mean(x,na.rm=T),1),"% ±",round(sd(x,na.rm=T),1),sep="") |
Also available in: Unified diff
updating figures