Revision d39ab57e
Added by Adam Wilson about 11 years ago
| climate/procedures/MOD35C5_Evaluation.r | ||
|---|---|---|
| 13 | 13 |
|
| 14 | 14 |
## get % cloudy |
| 15 | 15 |
mod09=raster("data/MOD09_2009.tif")
|
| 16 |
names(mod09)="MOD09CF" |
|
| 16 |
names(mod09)="C5MOD09CF"
|
|
| 17 | 17 |
NAvalue(mod09)=0 |
| 18 | 18 |
|
| 19 | 19 |
mod35c5=raster("data/MOD35_2009.tif")
|
| ... | ... | |
| 25 | 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 | 26 |
system("align.sh data/MOD35C6.vrt data/MOD09_2009.tif data/MOD35C6_2009.tif")
|
| 27 | 27 |
} |
| 28 |
mod35c6=raster("data/MOD35C6_2009.tif")
|
|
| 28 |
mod35c6=raster("data/MOD35C6_2009_v1.tif")
|
|
| 29 | 29 |
names(mod35c6)="C6MOD35CF" |
| 30 | 30 |
NAvalue(mod35c6)=255 |
| 31 | 31 |
|
| 32 | 32 |
## landcover |
| 33 | 33 |
if(!file.exists("data/MCD12Q1_IGBP_2009_051_wgs84_1km.tif")){
|
| 34 | 34 |
system(paste("/usr/local/gdal-1.10.0/bin/gdalwarp -tr 0.008983153 0.008983153 -r mode -ot Byte -co \"COMPRESS=LZW\"",
|
| 35 |
" /mnt/data/jetzlab/Data/environ/global/MODIS/MCD12Q1/051/MCD12Q1_051_2009.tif ", |
|
| 35 |
" /mnt/data/jetzlab/Data/environ/global/MODIS/MCD12Q1/051/MCD12Q1_051_2009_wgs84.tif ",
|
|
| 36 | 36 |
" -t_srs \"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs\" ", |
| 37 | 37 |
" -te -180.0044166 -60.0074610 180.0044166 90.0022083 ", |
| 38 | 38 |
"data/MCD12Q1_IGBP_2009_051_wgs84_1km.tif -overwrite ",sep=""))} |
| ... | ... | |
| 65 | 65 |
system("align.sh ~/acrobates/adamw/projects/interp/data/modis/MOD17/MOD17A3_Science_NPP_mean_00_12.tif data/MOD09_2009.tif data/MOD17.tif")
|
| 66 | 66 |
mod17=raster("data/MOD17.tif",format="GTiff")
|
| 67 | 67 |
NAvalue(mod17)=65535 |
| 68 |
names(mod17)="MOD17" |
|
| 68 |
names(mod17)="MOD17_unscaled"
|
|
| 69 | 69 |
|
| 70 | 70 |
if(!file.exists("data/MOD17qc.tif"))
|
| 71 | 71 |
system("align.sh ~/acrobates/adamw/projects/interp/data/modis/MOD17/MOD17A3_Science_NPP_Qc_mean_00_12.tif data/MOD09_2009.tif data/MOD17qc.tif")
|
| ... | ... | |
| 77 | 77 |
if(!file.exists("data/MOD11_2009.tif"))
|
| 78 | 78 |
system("align.sh ~/acrobates/adamw/projects/interp/data/modis/mod11/2009/MOD11_LST_2009.tif data/MOD09_2009.tif data/MOD11_2009.tif")
|
| 79 | 79 |
mod11=raster("data/MOD11_2009.tif",format="GTiff")
|
| 80 |
names(mod11)="MOD11" |
|
| 80 |
names(mod11)="MOD11_unscaled"
|
|
| 81 | 81 |
NAvalue(mod11)=0 |
| 82 | 82 |
if(!file.exists("data/MOD11qc_2009.tif"))
|
| 83 | 83 |
system("align.sh ~/acrobates/adamw/projects/interp/data/modis/mod11/2009/MOD11_Pmiss_2009.tif data/MOD09_2009.tif data/MOD11qc_2009.tif")
|
| 84 | 84 |
mod11qc=raster("data/MOD11qc_2009.tif",format="GTiff")
|
| 85 | 85 |
names(mod11qc)="MOD11CF" |
| 86 | 86 |
|
| 87 |
|
|
| 88 |
### Create some summary objects for plotting |
|
| 89 |
#difm=v6m-v5m |
|
| 90 |
#v5v6compare=stack(v5m,v6m,difm) |
|
| 91 |
#names(v5v6compare)=c("Collection 5","Collection 6","Difference (C6-C5)")
|
|
| 92 |
|
|
| 93 | 87 |
### Processing path |
| 94 | 88 |
if(!file.exists("data/MOD35pp.tif"))
|
| 95 | 89 |
system("align.sh data/MOD35_ProcessPath.tif data/MOD09_2009.tif data/MOD35pp.tif")
|
| ... | ... | |
| 115 | 109 |
bgyr=colorRampPalette(c("blue","green","yellow","red"))
|
| 116 | 110 |
cols=bgyr(n) |
| 117 | 111 |
|
| 118 |
#levelplot(lulcf,margin=F,layers="LULC") |
|
| 119 |
|
|
| 120 | 112 |
|
| 121 | 113 |
### Transects |
| 122 | 114 |
r1=Lines(list( |
| ... | ... | |
| 134 | 126 |
73.943,27.419, |
| 135 | 127 |
74.369,26.877 |
| 136 | 128 |
),ncol=2,byrow=T))),"India") |
| 137 |
r4=Lines(list( |
|
| 138 |
Line(matrix(c( |
|
| 139 |
-5.164,42.270, |
|
| 140 |
-4.948,42.162 |
|
| 141 |
),ncol=2,byrow=T))),"Spain") |
|
| 129 |
#r4=Lines(list(
|
|
| 130 |
# Line(matrix(c(
|
|
| 131 |
# -5.164,42.270,
|
|
| 132 |
# -4.948,42.162
|
|
| 133 |
# ),ncol=2,byrow=T))),"Spain")
|
|
| 142 | 134 |
|
| 143 | 135 |
r5=Lines(list( |
| 144 | 136 |
Line(matrix(c( |
| ... | ... | |
| 146 | 138 |
33.802,12.894 |
| 147 | 139 |
),ncol=2,byrow=T))),"Sudan") |
| 148 | 140 |
|
| 149 |
r6=Lines(list( |
|
| 150 |
Line(matrix(c( |
|
| 151 |
-63.353,-10.746, |
|
| 152 |
-63.376,-9.310 |
|
| 153 |
),ncol=2,byrow=T))),"Brazil") |
|
| 141 |
#r6=Lines(list(
|
|
| 142 |
# Line(matrix(c(
|
|
| 143 |
# -63.353,-10.746,
|
|
| 144 |
# -63.376,-9.310
|
|
| 145 |
# ),ncol=2,byrow=T))),"Brazil")
|
|
| 154 | 146 |
|
| 155 | 147 |
|
| 156 | 148 |
trans=SpatialLines(list(r1,r2,r3,r5),CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs "))
|
| ... | ... | |
| 179 | 171 |
td=crop(mod11,transb[x]) |
| 180 | 172 |
tdd=lapply(list(mod35c5,mod35c6,mod09,mod17,mod17qc,mod11,mod11qc,lulc,pp),function(l) resample(crop(l,transb[x]),td,method="ngb")) |
| 181 | 173 |
## normalize MOD11 and MOD17 |
| 182 |
for(j in which(do.call(c,lapply(tdd,function(i) names(i)))%in%c("MOD11","MOD17"))){
|
|
| 174 |
for(j in which(do.call(c,lapply(tdd,function(i) names(i)))%in%c("MOD11_unscaled","MOD17_unscaled"))){
|
|
| 183 | 175 |
trange=cellStats(tdd[[j]],range) |
| 184 | 176 |
tscaled=100*(tdd[[j]]-trange[1])/(trange[2]-trange[1]) |
| 185 | 177 |
tscaled@history=list(range=trange) |
| 186 |
names(tscaled)=paste(names(tdd[[j]]),"scaled",collapse="_",sep="_")
|
|
| 178 |
names(tscaled)=sub("_unscaled","",names(tdd[[j]]))
|
|
| 187 | 179 |
tdd=c(tdd,tscaled) |
| 188 | 180 |
} |
| 189 | 181 |
return(brick(tdd)) |
| ... | ... | |
| 218 | 210 |
|
| 219 | 211 |
|
| 220 | 212 |
## comparison of % cloudy days |
| 221 |
dif_c5_09=mod35c5-mod09 |
|
| 213 |
if(!file.exists("data/dif_c5_09.tif"))
|
|
| 214 |
overlay(mod35c5,mod09,fun=function(x,y) {return(x-y)},file="data/dif_c5_09.tif",format="GTiff",options=c("COMPRESS=LZW","ZLEVEL=9"),overwrite=T)
|
|
| 215 |
dif_c5_09=raster("data/dif_c5_09.tif",format="GTiff")
|
|
| 216 |
|
|
| 222 | 217 |
#dif_c6_09=mod35c6-mod09 |
| 223 | 218 |
#dif_c5_c6=mod35c5-mod35c6 |
| 224 | 219 |
|
| 225 |
## exploring various ways to compare cloud products |
|
| 226 |
t1=trd1[[1]] |
|
| 227 |
dif_p=calc(trd1[[1]], function(x) (x[1]-x[3])/(1-x[1])) |
|
| 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 |
|
|
| 232 |
names(edge)="edge" |
|
| 233 |
names(edgeb)="edgeb" |
|
| 234 |
|
|
| 235 |
td1=as.data.frame(stack(t1,edge,edgeb)) |
|
| 220 |
## exploring various ways to compare cloud products along landcover or processing path edges |
|
| 221 |
#t1=trd1[[1]] |
|
| 222 |
#dif_p=calc(trd1[[1]], function(x) (x[1]-x[3])/(1-x[1])) |
|
| 223 |
#edge=calc(edge(subset(t1,"MCD12Q1"),classes=T,type="inner"),function(x) ifelse(x==1,1,NA)) |
|
| 224 |
#edgeb=buffer(edge,width=5000) |
|
| 225 |
#edgeb=calc(edgeb,function(x) ifelse(is.na(x),0,1)) |
|
| 226 |
#names(edge)="edge" |
|
| 227 |
#names(edgeb)="edgeb" |
|
| 228 |
#td1=as.data.frame(stack(t1,edge,edgeb)) |
|
| 229 |
#cor(td1$MOD17,td1$C6MOD35,use="complete",method="spearman") |
|
| 230 |
#cor(td1$MOD17[td1$edgeb==1],td1$C5MOD35[td1$edgeb==1],use="complete",method="spearman") |
|
| 231 |
|
|
| 232 |
### Correlations |
|
| 233 |
#trdw=cast(trd,trans+x+y~variable,value="value") |
|
| 234 |
#cor(trdw$MOD17,trdw$C5MOD35,use="complete",method="spearman") |
|
| 236 | 235 |
|
| 237 |
cor(td1$MOD17,td1$C5MOD35,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") |
|
| 236 |
#Across all pixels in the four regions analyzed in Figure 3 there is a much larger correlation between mean NPP and the C5 MOD35 CF (Spearman’s ρ = -0.61, n=58,756) than the C6 MOD35 CF (ρ = 0.00, n=58,756) or MOD09 (ρ = -0.07, n=58,756) products. |
|
| 237 |
#by(trdw,trdw$trans,function(x) cor(as.data.frame(na.omit(x[,c("C5MOD35CF","C6MOD35CF","C5MOD09CF","MOD17","MOD11")])),use="complete",method="spearman"))
|
|
| 244 | 238 |
|
| 245 |
cor.test( |
|
| 246 |
splom(t1) |
|
| 247 |
plot(mod17,mod17qc) |
|
| 248 |
xyplot(MOD17~C5MOD35CF|edgeb,data=td1) |
|
| 249 |
bwplot(MCD12Q1~C5MOD35CF|edgeb,data=td1) |
|
| 250 | 239 |
|
| 251 |
plot(dif_p) |
|
| 240 |
## table of correlations |
|
| 241 |
#trdw_cor=as.data.frame(na.omit(trdw[,c("C5MOD35CF","C6MOD35CF","C5MOD09CF","MOD17","MOD11")]))
|
|
| 242 |
#nrow(trdw_cor) |
|
| 243 |
#round(cor(trdw_cor,method="spearman"),2) |
|
| 252 | 244 |
|
| 253 |
#rondonia=trd[trd$trans=="Brazil",] |
|
| 254 |
#trd=trd[trd$trans!="Brazil",] |
|
| 255 | 245 |
|
| 246 |
## set up some graphing parameters |
|
| 256 | 247 |
at=seq(0,100,leng=100) |
| 257 | 248 |
bgyr=colorRampPalette(c("purple","blue","green","yellow","orange","red","red"))
|
| 258 | 249 |
bgrayr=colorRampPalette(c("purple","blue","grey","red","red"))
|
| ... | ... | |
| 265 | 256 |
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 | 257 |
|
| 267 | 258 |
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))
|
| 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)) |
|
| 259 |
g2$strip=strip.custom(var.name="Difference (C5MOD35-C5MOD09)",style=1,strip.names=T,strip.levels=F) #update strip text |
|
| 260 |
#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)) |
|
| 273 | 261 |
|
| 274 | 262 |
### regional plots |
| 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"),
|
|
| 263 |
p1=useOuterStrips(levelplot(value~x*y|variable+trans,data=trd[!trd$variable%in%c("MOD17_unscaled","MOD11_unscaled","MCD12Q1","MOD35pp"),],asp=1,scales=list(draw=F,rot=0,relation="free"),
|
|
| 276 | 264 |
at=at,col.regions=cols,maxpixels=7e6, |
| 277 |
ylab="Latitude",xlab="Longitude"),strip = strip.custom(par.strip.text=list(cex=.75)))+layer(sp.lines(trans,lwd=2))
|
|
| 265 |
ylab="Latitude",xlab="Longitude"),strip = strip.custom(par.strip.text=list(cex=.7)))+layer(sp.lines(trans,lwd=2)) |
|
| 278 | 266 |
|
| 279 | 267 |
p2=useOuterStrips( |
| 280 | 268 |
levelplot(value~x*y|variable+trans,data=trd[trd$variable%in%c("MCD12Q1"),],
|
| ... | ... | |
| 284 | 272 |
right=list(fun=draw.key(list(columns=1,#title="MCD12Q1 \n IGBP Land \n Cover", |
| 285 | 273 |
rectangles=list(col=IGBP$col,size=1), |
| 286 | 274 |
text=list(as.character(IGBP$ID),at=IGBP$ID-.5))))), |
| 287 |
ylab="",xlab=" "),strip = strip.custom(par.strip.text=list(cex=.75)),strip.left=F)+layer(sp.lines(trans,lwd=2))
|
|
| 275 |
ylab="",xlab=" "),strip = strip.custom(par.strip.text=list(cex=.7)),strip.left=F)+layer(sp.lines(trans,lwd=2)) |
|
| 288 | 276 |
p3=useOuterStrips( |
| 289 | 277 |
levelplot(value~x*y|variable+trans,data=trd[trd$variable%in%c("MOD35pp"),],
|
| 290 | 278 |
asp=1,scales=list(draw=F,rot=0,relation="free"),colorkey=F, |
| ... | ... | |
| 293 | 281 |
right=list(fun=draw.key(list(columns=1,#title="MOD35 \n Processing \n Path", |
| 294 | 282 |
rectangles=list(col=c("blue","cyan","tan","darkgreen"),size=1),
|
| 295 | 283 |
text=list(c("Water","Coast","Desert","Land")))))),
|
| 296 |
ylab="",xlab=" "),strip = strip.custom(par.strip.text=list(cex=.75)),strip.left=F)+layer(sp.lines(trans,lwd=2))
|
|
| 284 |
ylab="",xlab=" "),strip = strip.custom(par.strip.text=list(cex=.7)),strip.left=F)+layer(sp.lines(trans,lwd=2)) |
|
| 297 | 285 |
|
| 298 | 286 |
## transects |
| 299 | 287 |
p4=xyplot(value~dist|transect,groups=variable,type=c("smooth","p"),
|
| 300 | 288 |
data=transd,panel=function(...,subscripts=subscripts) {
|
| 301 | 289 |
td=transd[subscripts,] |
| 302 | 290 |
## mod09 |
| 303 |
imod09=td$variable=="MOD09CF" |
|
| 291 |
imod09=td$variable=="C5MOD09CF"
|
|
| 304 | 292 |
panel.xyplot(td$dist[imod09],td$value[imod09],type=c("p","smooth"),span=0.2,subscripts=1:sum(imod09),col="red",pch=16,cex=.25)
|
| 305 | 293 |
## mod35C5 |
| 306 | 294 |
imod35=td$variable=="C5MOD35CF" |
| ... | ... | |
| 323 | 311 |
panel.xyplot(td$dist[imod11qc],td$value[imod11qc],type=c("p","smooth"),npoints=100,span=qcspan,subscripts=1:sum(imod11qc),col="orange",pch=16,cex=.25)
|
| 324 | 312 |
## land |
| 325 | 313 |
path=td[td$variable=="MOD35pp",] |
| 326 |
panel.segments(path$dist,-5,c(path$dist[-1],max(path$dist,na.rm=T)),-5,col=c("blue","cyan","tan","darkgreen")[path$value+1],subscripts=1:nrow(path),lwd=15,type="l")
|
|
| 314 |
panel.segments(path$dist,-10,c(path$dist[-1],max(path$dist,na.rm=T)),-10,col=c("blue","cyan","tan","darkgreen")[path$value+1],subscripts=1:nrow(path),lwd=10,type="l")
|
|
| 327 | 315 |
land=td[td$variable=="MCD12Q1",] |
| 328 |
panel.segments(land$dist,-10,c(land$dist[-1],max(land$dist,na.rm=T)),-10,col=IGBP$col[land$value+1],subscripts=1:nrow(land),lwd=15,type="l")
|
|
| 316 |
panel.segments(land$dist,-20,c(land$dist[-1],max(land$dist,na.rm=T)),-20,col=IGBP$col[land$value+1],subscripts=1:nrow(land),lwd=10,type="l")
|
|
| 329 | 317 |
},subscripts=T,par.settings = list(grid.pars = list(lineend = "butt")), |
| 330 | 318 |
scales=list( |
| 331 | 319 |
x=list(alternating=1,relation="free"),#, lim=c(0,70)), |
| 332 |
y=list(at=c(-10,-5,seq(0,100,len=5)), |
|
| 333 |
labels=c("IGBP","MOD35",seq(0,100,len=5)),
|
|
| 334 |
lim=c(-15,100))), |
|
| 320 |
y=list(at=c(-18,-10,seq(0,100,len=5)), |
|
| 321 |
labels=c("MCD12Q1 IGBP","MOD35 path",seq(0,100,len=5)),
|
|
| 322 |
lim=c(-25,100)), |
|
| 323 |
alternating=F), |
|
| 335 | 324 |
xlab="Distance Along Transect (km)", ylab="% Missing Data / % of Maximum Value", |
| 336 | 325 |
legend=list( |
| 337 | 326 |
bottom=list(fun=draw.key(list( rep=FALSE,columns=1,title=" ", |
| 338 |
## text=list(c("MOD09 % Cloudy","C5 MOD35 % Cloudy","C6 MOD35 % Cloudy","MOD17 % Missing","MOD17 (scaled)","MOD11 % Missing","MOD11 (scaled)")),
|
|
| 339 |
lines=list(type=c("b","b","b","b","b","l","b","l"),pch=16,cex=.5,lty=c(0,1,1,1,1,5,1,5),col=c("transparent","red","blue","black","darkgreen","darkgreen","orange","orange")),
|
|
| 340 |
text=list(c("MODIS Products","MOD09 % Cloudy","C5 MOD35 % Cloudy","C6 MOD35 % Cloudy","MOD17 % Missing","MOD17 (scaled)","MOD11 % Missing","MOD11 (scaled)")),
|
|
| 327 |
lines=list(type=c("b","b","b","b","b","l","b","l"),pch=16,cex=.5,
|
|
| 328 |
lty=c(0,1,1,1,1,5,1,5), |
|
| 329 |
col=c("transparent","red","blue","black","darkgreen","darkgreen","orange","orange")),
|
|
| 330 |
text=list( |
|
| 331 |
c("MODIS Products","C5 MOD09 % Cloudy","C5 MOD35 % Cloudy","C6 MOD35 % Cloudy","MOD17 % Missing","MOD17 (scaled)","MOD11 % Missing","MOD11 (scaled)")),
|
|
| 341 | 332 |
rectangles=list(border=NA,col=c(NA,"tan","darkgreen")), |
| 342 | 333 |
text=list(c("C5 MOD35 Processing Path","Desert","Land")),
|
| 343 |
rectangles=list(border=NA,col=c(NA,IGBP$col[sort(unique(transd$value[transd$variable=="MCD12Q1"]+1))])), |
|
| 344 |
text=list(c("MCD12Q1 IGBP Land Cover",IGBP$class[sort(unique(transd$value[transd$variable=="MCD12Q1"]+1))])))))),
|
|
| 345 |
strip = strip.custom(par.strip.text=list(cex=.75))) |
|
| 346 |
#print(p4) |
|
| 347 |
|
|
| 348 |
#trdw=cast(trd,trans+x+y~variable,value="value") |
|
| 349 |
#transdw=cast(transd,transect+dist~variable,value="value") |
|
| 350 |
#xyplot(MOD11CF~C5MOD35CF|transect,groups=MCD12Q1,data=transdw,pch=16,cex=1,scales=list(relation="free")) |
|
| 351 |
#xyplot(MOD17~C5MOD35CF|trans,groups=MCD12Q1,data=trdw,pch=16,cex=1,scales=list(relation="free")) |
|
| 334 |
rectangles=list(border=NA,col=c(NA,IGBP$col[sort(unique(transd$value[transd$variable=="MCD12Q1"]+1))])), |
|
| 335 |
text=list(c("MCD12Q1 IGBP Land Cover",IGBP$class[sort(unique(transd$value[transd$variable=="MCD12Q1"]+1))])))))),
|
|
| 336 |
strip = strip.custom(par.strip.text=list(cex=.75))) |
|
| 337 |
print(p4) |
|
| 352 | 338 |
|
| 353 |
#p5=levelplot(value~x*y|variable,data=rondonia,asp=1,scales=list(draw=F,rot=0,relation="free"),colorkey=T)#, |
|
| 354 |
#print(p5) |
|
| 355 | 339 |
|
| 356 | 340 |
|
| 357 | 341 |
CairoPDF("output/mod35compare.pdf",width=11,height=7)
|
| 358 | 342 |
#CairoPNG("output/mod35compare_%d.png",units="in", width=11,height=8.5,pointsize=4000,dpi=1200,antialias="subpixel")
|
| 359 | 343 |
### Global Comparison |
| 360 |
print(g1) |
|
| 361 | 344 |
print(g1,position=c(0,.35,1,1),more=T) |
| 362 | 345 |
print(g2,position=c(0,0,1,0.415),more=F) |
| 363 | 346 |
#print(g3,position=c(0.31,0.06,.42,0.27),more=F) |
| ... | ... | |
| 393 | 376 |
td[td$transect=="Venezuela",] |
| 394 | 377 |
|
| 395 | 378 |
|
| 396 |
|
|
| 397 |
|
|
| 398 |
## scatterplot of MOD35 vs MOD09 |
|
| 399 |
trdl=cast(trans+x+y~variable,value="value",data=trd) |
|
| 400 |
xyplot(MOD35C5qc~MOD09qc|trans+as.factor(MOD35pp),pch=16,cex=.2,data=trdl,auto.key=T)+layer(panel.abline(0,1)) |
|
| 401 |
|
|
| 402 |
|
|
| 403 |
### LANDCOVER |
|
| 404 |
levelplot(lulc,col.regions=IGBP$col,scales=list(cex=2),colorkey=list(space="right",at=0:17,labels=list(at=seq(0.5,16.5,by=1),labels=levels(lulc)[[1]]$class,cex=2)),margin=F) |
|
| 405 |
|
|
| 406 |
levelplot(mcompare,col.regions=cols,at=at,margin=F,sub="Frequency of MOD35 Clouds in March") |
|
| 407 |
#levelplot(dif,col.regions=bgyr(20),margin=F) |
|
| 408 |
levelplot(mdiff,col.regions=bgyr(100),at=seq(mdiff@data@min,mdiff@data@max,len=100),margin=F) |
|
| 409 |
|
|
| 410 |
|
|
| 411 |
boxplot(as.matrix(subset(dif,subset=1))~forest,varwidth=T,notch=T);abline(h=0) |
|
| 412 |
|
|
| 413 |
|
|
| 414 |
levelplot(modprod,main="Missing Data (%) in MOD17 (NPP) and MOD43 (BRDF Reflectance)", |
|
| 415 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at) |
|
| 416 |
|
|
| 417 |
|
|
| 418 |
|
|
| 419 |
|
|
| 420 |
levelplot(modprod,main="Missing Data (%) in MOD17 (NPP) and MOD43 (BRDF Reflectance)", |
|
| 421 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at, |
|
| 422 |
xlim=c(-7300000,-6670000),ylim=c(0,600000)) |
|
| 423 |
|
|
| 424 |
levelplot(v5m,main="Missing Data (%) in MOD17 (NPP) and MOD43 (BRDF Reflectance)", |
|
| 425 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at, |
|
| 426 |
xlim=c(-7200000,-6670000),ylim=c(0,400000),margin=F) |
|
| 427 |
|
|
| 428 |
|
|
| 429 |
levelplot(subset(v5v6compare,1:2),main="Proportion Cloudy Days (%) in Collection 5 and 6 MOD35", |
|
| 430 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at, |
|
| 431 |
margin=F) |
|
| 432 |
|
|
| 433 |
levelplot(subset(v5v6compare,1:2),main="Proportion Cloudy Days (%) in Collection 5 and 6 MOD35", |
|
| 434 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at, |
|
| 435 |
xlim=c(-7200000,-6670000),ylim=c(0,400000),margin=F) |
|
| 436 |
|
|
| 437 |
levelplot(subset(v5v6compare,1:2),main="Proportion Cloudy Days (%) in Collection 5 and 6 MOD35", |
|
| 438 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at, |
|
| 439 |
xlim=c(-7500000,-7200000),ylim=c(700000,1000000),margin=F) |
|
| 440 |
|
|
| 441 |
|
|
| 442 |
dev.off() |
|
| 443 |
|
|
| 444 |
### smoothing plots |
|
| 445 |
## explore smoothed version |
|
| 446 |
td=subset(v6,m) |
|
| 447 |
## build weight matrix |
|
| 448 |
s=3 |
|
| 449 |
w=matrix(1/(s*s),nrow=s,ncol=s) |
|
| 450 |
#w[s-1,s-1]=4/12; w |
|
| 451 |
td2=focal(td,w=w) |
|
| 452 |
td3=stack(td,td2) |
|
| 453 |
|
|
| 454 |
levelplot(td3,col.regions=cols,at=at,margin=F) |
|
| 455 |
|
|
| 456 |
dev.off() |
|
| 457 |
plot(stack(difm,lulc)) |
|
| 458 |
|
|
| 459 |
### ROI |
|
| 460 |
tile_ll=projectExtent(v6, "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") |
|
| 461 |
|
|
| 462 |
62,59 |
|
| 463 |
0,3 |
|
| 464 |
|
|
| 465 |
|
|
| 466 |
|
|
| 467 |
#### export KML timeseries |
|
| 379 |
#### export KML |
|
| 468 | 380 |
library(plotKML) |
| 469 |
tile="h11v08" |
|
| 470 |
file=paste("summary/MOD35_",tile,".nc",sep="")
|
|
| 471 |
system(paste("gdalwarp -overwrite -multi -ot INT16 -r cubicspline -srcnodata 255 -dstnodata 255 -s_srs '+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs' -t_srs 'EPSG:4326' NETCDF:",file,":PCloud MOD35_",tile,".tif",sep=""))
|
|
| 472 |
|
|
| 473 |
v6sp=brick(paste("MOD35_",tile,".tif",sep=""))
|
|
| 474 |
v6sp=readAll(v6sp) |
|
| 475 | 381 |
|
| 476 |
## wasn't working with line below, perhaps Z should just be text? not date? |
|
| 477 |
v6sp=setZ(v6sp,as.Date(paste("2011-",1:12,"-15",sep="")))
|
|
| 478 |
names(v6sp)=month.name |
|
| 382 |
kml_open("output/modiscloud.kml")
|
|
| 479 | 383 |
|
| 480 |
kml_open("output/mod35.kml")
|
|
| 384 |
readAll(mod35c5)
|
|
| 481 | 385 |
|
| 482 |
|
|
| 483 |
kml_layer.RasterBrick(v6sp, |
|
| 484 |
plot.legend = TRUE, dtime = "", tz = "GMT", |
|
| 485 |
z.lim = c(0,100),colour_scale = get("colour_scale_numeric", envir = plotKML.opts))
|
|
| 386 |
kml_layer.Raster(mod35c5, |
|
| 387 |
plot.legend = TRUE,raster_name="Collection 5 MOD35 Cloud Frequency", |
|
| 388 |
z.lim = c(0,100),colour_scale = get("colour_scale_numeric", envir = plotKML.opts),
|
|
| 486 | 389 |
# home_url = get("home_url", envir = plotKML.opts),
|
| 487 | 390 |
# metadata = NULL, html.table = NULL, |
| 488 |
# altitudeMode = "clampToGround", balloon = FALSE,
|
|
| 391 |
altitudeMode = "clampToGround", balloon = FALSE
|
|
| 489 | 392 |
) |
| 490 | 393 |
|
| 394 |
system(paste("gdal_translate -of KMLSUPEROVERLAY ",mod35c5@file@name," output/mod35c5.kmz -co FORMAT=JPEG"))
|
|
| 395 |
|
|
| 491 | 396 |
logo = "http://static.tumblr.com/t0afs9f/KWTm94tpm/yale_logo.png" |
| 492 | 397 |
kml_screen(image.file = logo, position = "UL", sname = "YALE logo",size=c(.1,.1)) |
| 493 |
kml_close("mod35.kml")
|
|
| 494 |
kml_compress("mod35.kml",files=c(paste(month.name,".png",sep=""),"obj_legend.png"),zip="/usr/bin/zip")
|
|
| 398 |
kml_close("modiscloud.kml")
|
|
| 399 |
kml_compress("modiscloud.kml",files=c(paste(month.name,".png",sep=""),"obj_legend.png"),zip="/usr/bin/zip")
|
|
Also available in: Unified diff
Submitted MOD35-landcover bias paper with code from this commit. Also added short script to test swtif program.