Revision 14ad4a96
Added by Adam Wilson over 11 years ago
| climate/procedures/MOD35C5_Evaluation.r | ||
|---|---|---|
| 9 | 9 |
library(Cairo) |
| 10 | 10 |
library(reshape) |
| 11 | 11 |
library(rgeos) |
| 12 |
library(splancs) |
|
| 12 | 13 |
|
| 13 | 14 |
## get % cloudy |
| 14 | 15 |
mod09=raster("data/MOD09_2009.tif")
|
| ... | ... | |
| 21 | 22 |
|
| 22 | 23 |
## mod35C6 annual |
| 23 | 24 |
if(!file.exists("data/MOD35C6_2009.tif")){
|
| 24 |
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/*2009mean.nc ")
|
|
| 25 |
system("align.sh data/MOD35C6.vrt data/MOD09_2009.tif 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")
|
|
| 26 | 27 |
} |
| 27 |
mod35c6=raster("data/MOD35C6_2009.tif")
|
|
| 28 |
mod35c6=raster("data/MOD35C6_2009_v1.tif")
|
|
| 28 | 29 |
names(mod35c6)="C6MOD35CF" |
| 29 | 30 |
NAvalue(mod35c6)=255 |
| 30 | 31 |
|
| ... | ... | |
| 158 | 159 |
writeOGR(SpatialLinesDataFrame(trans,data=data.frame(ID=names(trans)),match.ID=F),"output",layer="transects",driver="ESRI Shapefile",overwrite=T) |
| 159 | 160 |
|
| 160 | 161 |
## buffer transects to get regional values |
| 161 |
transb=gBuffer(trans,0.4) |
|
| 162 |
transb=gBuffer(trans,byid=T,width=0.4)
|
|
| 162 | 163 |
|
| 163 | 164 |
## make polygons of bounding boxes |
| 164 | 165 |
bb0 <- lapply(slot(transb, "polygons"), bbox) |
| 165 |
library(splancs) |
|
| 166 | 166 |
bb1 <- lapply(bb0, bboxx) |
| 167 | 167 |
# turn these into matrices using a helper function in splancs |
| 168 | 168 |
bb2 <- lapply(bb1, function(x) rbind(x, x[1,])) |
| ... | ... | |
| 176 | 176 |
# loop over the closed matrix rings, adding the IDs |
| 177 | 177 |
bbs <- SpatialPolygons(bb3, proj4string=CRS(proj4string(transb))) |
| 178 | 178 |
|
| 179 |
|
|
| 180 | 179 |
trd1=lapply(1:length(transb),function(x) {
|
| 181 | 180 |
td=crop(mod11,transb[x]) |
| 182 | 181 |
tdd=lapply(list(mod35c5,mod35c6,mod09,mod17,mod17qc,mod11,mod11qc,lulc,pp),function(l) resample(crop(l,transb[x]),td,method="ngb")) |
| ... | ... | |
| 219 | 218 |
transd$type=ifelse(grepl("MOD35|MOD09|CF",transd$variable),"CF","Data")
|
| 220 | 219 |
|
| 221 | 220 |
|
| 222 |
|
|
| 223 | 221 |
## comparison of % cloudy days |
| 224 | 222 |
dif_c5_09=mod35c5-mod09 |
| 225 |
dif_c6_09=mod35c6-mod09 |
|
| 226 |
dif_c5_c6=mod35c5-mod35c6 |
|
| 223 |
#dif_c6_09=mod35c6-mod09
|
|
| 224 |
#dif_c5_c6=mod35c5-mod35c6
|
|
| 227 | 225 |
|
| 226 |
## exploring various ways to compare cloud products |
|
| 228 | 227 |
t1=trd1[[1]] |
| 229 | 228 |
dif_p=calc(trd1[[1]], function(x) (x[1]-x[3])/(1-x[1])) |
| 230 | 229 |
edge=edge(subset(t1,"MCD12Q1"),classes=T,type="inner") |
| 231 | 230 |
names(edge)="edge" |
| 232 | 231 |
td1=as.data.frame(stack(t1,edge)) |
| 233 | 232 |
|
| 234 |
|
|
| 235 | 233 |
cor(td1$MOD17,td1$C5MOD35,use="complete",method="spearman") |
| 236 | 234 |
cor(td1$MOD17[td1$edge==1],td1$C5MOD35[td1$edge==1],use="complete",method="spearman") |
| 237 |
|
|
| 238 | 235 |
cor(td1,use="complete",method="spearman") |
| 239 |
|
|
| 240 | 236 |
splom(t1) |
| 241 |
|
|
| 242 | 237 |
plot(mod17,mod17qc) |
| 243 |
|
|
| 244 | 238 |
xyplot(MOD17~C5MOD35CF|edge,data=td1) |
| 245 |
|
|
| 246 |
, function(x) (x[1]-x[3])/(1-x[1])) |
|
| 247 |
|
|
| 248 | 239 |
plot(dif_p) |
| 249 | 240 |
|
| 250 | 241 |
#rondonia=trd[trd$trans=="Brazil",] |
| ... | ... | |
| 260 | 251 |
coast=map2SpatialLines(map("world", interior=FALSE, plot=FALSE),proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
|
| 261 | 252 |
|
| 262 | 253 |
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)) |
| 263 |
g2=levelplot(dif,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))
|
|
| 254 |
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))
|
|
| 264 | 255 |
trellis.par.set(background=list(fill="white"),panel.background=list(fill="white")) |
| 265 |
#g3=histogram(dif,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))
|
|
| 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))
|
|
| 266 | 257 |
|
| 267 | 258 |
### regional plots |
| 268 | 259 |
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"),
|
| ... | ... | |
| 338 | 329 |
strip = strip.custom(par.strip.text=list(cex=.75))) |
| 339 | 330 |
print(p4) |
| 340 | 331 |
|
| 341 |
trdw=cast(trd,trans+x+y~variable,value="value") |
|
| 342 |
transdw=cast(transd,transect+dist~variable,value="value") |
|
| 343 |
|
|
| 344 |
xyplot(MOD11CF~C5MOD35CF|transect,groups=MCD12Q1,data=transdw,pch=16,cex=1,scales=list(relation="free")) |
|
| 345 |
xyplot(MOD17~C5MOD35CF|trans,groups=MCD12Q1,data=trdw,pch=16,cex=1,scales=list(relation="free")) |
|
| 332 |
#trdw=cast(trd,trans+x+y~variable,value="value") |
|
| 333 |
#transdw=cast(transd,transect+dist~variable,value="value") |
|
| 334 |
#xyplot(MOD11CF~C5MOD35CF|transect,groups=MCD12Q1,data=transdw,pch=16,cex=1,scales=list(relation="free")) |
|
| 335 |
#xyplot(MOD17~C5MOD35CF|trans,groups=MCD12Q1,data=trdw,pch=16,cex=1,scales=list(relation="free")) |
|
| 346 | 336 |
|
| 347 | 337 |
#p5=levelplot(value~x*y|variable,data=rondonia,asp=1,scales=list(draw=F,rot=0,relation="free"),colorkey=T)#, |
| 348 | 338 |
#print(p5) |
| ... | ... | |
| 352 | 342 |
#CairoPNG("output/mod35compare_%d.png",units="in", width=11,height=8.5,pointsize=4000,dpi=1200,antialias="subpixel")
|
| 353 | 343 |
### Global Comparison |
| 354 | 344 |
print(g1) |
| 355 |
#print(g1,position=c(0,.33,1,1),more=T)
|
|
| 356 |
#print(g2,position=c(0,0,1,0.394),more=T)
|
|
| 357 |
#print(g3,position=c(0.31,0.06,.42,0.27),more=F)
|
|
| 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) |
|
| 358 | 348 |
### MOD35 Desert Processing path |
| 359 | 349 |
levelplot(pp,asp=1,scales=list(draw=T,rot=0),maxpixels=1e6, |
| 360 | 350 |
at=c(-1:3),col.regions=c("blue","cyan","tan","darkgreen"),margin=F,
|
Also available in: Unified diff
Updating figures for C5MOD35 Manuscript