Revision 5f212fe8
Added by Adam Wilson over 11 years ago
| climate/procedures/MOD35C5_Evaluation.r | ||
|---|---|---|
| 25 | 25 |
system("align.sh data/MOD35C6.vrt data/MOD09_2009.tif data/MOD35C6_2009.tif")
|
| 26 | 26 |
} |
| 27 | 27 |
mod35c6=raster("data/MOD35C6_2009.tif")
|
| 28 |
mod35c6=crop(mod35c6,mod09) |
|
| 29 | 28 |
names(mod35c6)="C6MOD35CF" |
| 30 | 29 |
NAvalue(mod35c6)=255 |
| 31 | 30 |
|
| 32 | 31 |
|
| 33 | 32 |
## landcover |
| 34 |
if(!file.exists("data/MCD12Q1_IGBP_2009_v5_wgs84_1km.tif")){
|
|
| 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\"",
|
| 36 |
" ~/acrobatesroot/jetzlab/Data/environ/global/landcover/MODIS/MCD12Q1_IGBP_2009_v5_wgs84.tif ",
|
|
| 35 |
" /mnt/data/jetzlab/Data/environ/global/MODIS/MCD12Q1/051/MCD12Q1_051_2009.tif ",
|
|
| 37 | 36 |
" -t_srs \"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs\" ", |
| 38 |
"data/MCD12Q1_IGBP_2009_v5_wgs84_1km.tif -overwrite ",sep=""))} |
|
| 39 |
lulc=raster("data/MCD12Q1_IGBP_2009_v5_wgs84_1km.tif")
|
|
| 37 |
" -te -180.0044166 -60.0074610 180.0044166 90.0022083 ", |
|
| 38 |
"data/MCD12Q1_IGBP_2009_051_wgs84_1km.tif -overwrite ",sep=""))} |
|
| 39 |
lulc=raster("data/MCD12Q1_IGBP_2009_051_wgs84_1km.tif")
|
|
| 40 | 40 |
|
| 41 |
## read it in |
|
| 42 |
lulc=raster("data/MCD12Q1_IGBP_2005_v51_1km_wgs84.tif")
|
|
| 43 | 41 |
# lulc=ratify(lulc) |
| 44 | 42 |
data(worldgrids_pal) #load palette |
| 45 | 43 |
IGBP=data.frame(ID=0:16,col=worldgrids_pal$IGBP[-c(18,19)], |
| 46 | 44 |
lulc_levels2=c("Water","Forest","Forest","Forest","Forest","Forest","Shrublands","Shrublands","Savannas","Savannas","Grasslands","Permanent wetlands","Croplands","Urban and built-up","Cropland/Natural vegetation mosaic","Snow and ice","Barren or sparsely vegetated"),stringsAsFactors=F)
|
| 47 | 45 |
IGBP$class=rownames(IGBP);rownames(IGBP)=1:nrow(IGBP) |
| 48 | 46 |
levels(lulc)=list(IGBP) |
| 49 |
lulc=crop(lulc,mod09) |
|
| 47 |
#lulc=crop(lulc,mod09)
|
|
| 50 | 48 |
names(lulc)="MCD12Q1" |
| 51 | 49 |
|
| 52 | 50 |
## make land mask |
| ... | ... | |
| 54 | 52 |
land=calc(lulc,function(x) ifelse(x==0,NA,1),file="data/land.tif",options=c("COMPRESS=LZW","ZLEVEL=9","PREDICTOR=2"),datatype="INT1U",overwrite=T)
|
| 55 | 53 |
land=raster("data/land.tif")
|
| 56 | 54 |
|
| 57 |
|
|
| 58 | 55 |
## mask cloud masks to land pixels |
| 59 | 56 |
#mod09l=mask(mod09,land) |
| 60 | 57 |
#mod35l=mask(mod35,land) |
| ... | ... | |
| 68 | 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")
|
| 69 | 66 |
mod17=raster("data/MOD17.tif",format="GTiff")
|
| 70 | 67 |
NAvalue(mod17)=65535 |
| 71 |
mod17=crop(mod17,mod09) |
|
| 72 | 68 |
names(mod17)="MOD17" |
| 73 | 69 |
|
| 74 | 70 |
if(!file.exists("data/MOD17qc.tif"))
|
| 75 | 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")
|
| 76 | 72 |
mod17qc=raster("data/MOD17qc.tif",format="GTiff")
|
| 77 | 73 |
NAvalue(mod17qc)=255 |
| 78 |
mod17qc=crop(mod17qc,mod09) |
|
| 79 | 74 |
names(mod17qc)="MOD17CF" |
| 80 | 75 |
|
| 81 | 76 |
## MOD11 via earth engine |
| ... | ... | |
| 87 | 82 |
if(!file.exists("data/MOD11qc_2009.tif"))
|
| 88 | 83 |
system("align.sh ~/acrobates/adamw/projects/interp/data/modis/mod11/2009/MOD11_Pmiss_2009.tif data/MOD09_2009.tif data/MOD11qc_2009.tif")
|
| 89 | 84 |
mod11qc=raster("data/MOD11qc_2009.tif",format="GTiff")
|
| 90 |
mod11qc=crop(mod11qc,mod09) |
|
| 91 | 85 |
names(mod11qc)="MOD11CF" |
| 92 | 86 |
|
| 93 | 87 |
|
| ... | ... | |
| 102 | 96 |
pp=raster("data/MOD35pp.tif")
|
| 103 | 97 |
NAvalue(pp)=255 |
| 104 | 98 |
names(pp)="MOD35pp" |
| 105 |
pp=crop(pp,mod09) |
|
| 106 | 99 |
|
| 107 |
## comparison of % cloudy days |
|
| 108 |
dif_c5_09=mod35c5-mod09 |
|
| 109 |
dif_c6_09=mod35c6-mod09 |
|
| 110 |
dif_c5_c6=mod35c5-mod35c6 |
|
| 111 | 100 |
|
| 112 | 101 |
#hist(dif,maxsamp=1000000) |
| 113 | 102 |
## draw lulc-stratified random sample of mod35-mod09 differences |
| ... | ... | |
| 194 | 183 |
## normalize MOD11 and MOD17 |
| 195 | 184 |
for(j in which(do.call(c,lapply(tdd,function(i) names(i)))%in%c("MOD11","MOD17"))){
|
| 196 | 185 |
trange=cellStats(tdd[[j]],range) |
| 197 |
tdd[[j]]=100*(tdd[[j]]-trange[1])/(trange[2]-trange[1]) |
|
| 198 |
tdd[[j]]@history=list(range=trange) |
|
| 186 |
tscaled=100*(tdd[[j]]-trange[1])/(trange[2]-trange[1]) |
|
| 187 |
tscaled@history=list(range=trange) |
|
| 188 |
names(tscaled)=paste(names(tdd[[j]]),"scaled",collapse="_",sep="_") |
|
| 189 |
tdd=c(tdd,tscaled) |
|
| 199 | 190 |
} |
| 200 | 191 |
return(brick(tdd)) |
| 201 | 192 |
}) |
| ... | ... | |
| 227 | 218 |
|
| 228 | 219 |
transd$type=ifelse(grepl("MOD35|MOD09|CF",transd$variable),"CF","Data")
|
| 229 | 220 |
|
| 221 |
|
|
| 222 |
|
|
| 223 |
## comparison of % cloudy days |
|
| 224 |
dif_c5_09=mod35c5-mod09 |
|
| 225 |
dif_c6_09=mod35c6-mod09 |
|
| 226 |
dif_c5_c6=mod35c5-mod35c6 |
|
| 227 |
|
|
| 228 |
t1=trd1[[1]] |
|
| 229 |
dif_p=calc(trd1[[1]], function(x) (x[1]-x[3])/(1-x[1])) |
|
| 230 |
edge=edge(subset(t1,"MCD12Q1"),classes=T,type="inner") |
|
| 231 |
names(edge)="edge" |
|
| 232 |
td1=as.data.frame(stack(t1,edge)) |
|
| 233 |
|
|
| 234 |
|
|
| 235 |
cor(td1$MOD17,td1$C5MOD35,use="complete",method="spearman") |
|
| 236 |
cor(td1$MOD17[td1$edge==1],td1$C5MOD35[td1$edge==1],use="complete",method="spearman") |
|
| 237 |
|
|
| 238 |
cor(td1,use="complete",method="spearman") |
|
| 239 |
|
|
| 240 |
splom(t1) |
|
| 241 |
|
|
| 242 |
plot(mod17,mod17qc) |
|
| 243 |
|
|
| 244 |
xyplot(MOD17~C5MOD35CF|edge,data=td1) |
|
| 245 |
|
|
| 246 |
, function(x) (x[1]-x[3])/(1-x[1])) |
|
| 247 |
|
|
| 248 |
plot(dif_p) |
|
| 249 |
|
|
| 230 | 250 |
#rondonia=trd[trd$trans=="Brazil",] |
| 231 | 251 |
#trd=trd[trd$trans!="Brazil",] |
| 232 | 252 |
|
| ... | ... | |
| 239 | 259 |
library(maptools) |
| 240 | 260 |
coast=map2SpatialLines(map("world", interior=FALSE, plot=FALSE),proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
|
| 241 | 261 |
|
| 242 |
g1=levelplot(stack(mod35c5,mod35c6,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))
|
|
| 243 |
#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))
|
|
| 244 |
#trellis.par.set(background=list(fill="white"),panel.background=list(fill="white"))
|
|
| 262 |
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))
|
|
| 264 |
trellis.par.set(background=list(fill="white"),panel.background=list(fill="white")) |
|
| 245 | 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)) |
| 246 | 266 |
|
| 247 | 267 |
### regional plots |
| ... | ... | |
| 309 | 329 |
legend=list( |
| 310 | 330 |
bottom=list(fun=draw.key(list( rep=FALSE,columns=1,title=" ", |
| 311 | 331 |
## text=list(c("MOD09 % Cloudy","C5 MOD35 % Cloudy","C6 MOD35 % Cloudy","MOD17 % Missing","MOD17 (scaled)","MOD11 % Missing","MOD11 (scaled)")),
|
| 312 |
lines=list(type=c("b","b","b","b","l","b","l"),pch=16,cex=.5,lty=c(1,1,1,1,5,1,5),col=c("red","blue","black","darkgreen","darkgreen","orange","orange")),
|
|
| 313 |
text=list(c("MOD09 % Cloudy","C5 MOD35 % Cloudy","C6 MOD35 % Cloudy","MOD17 % Missing","MOD17 (scaled)","MOD11 % Missing","MOD11 (scaled)")),
|
|
| 332 |
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")),
|
|
| 333 |
text=list(c("MODIS Products","MOD09 % Cloudy","C5 MOD35 % Cloudy","C6 MOD35 % Cloudy","MOD17 % Missing","MOD17 (scaled)","MOD11 % Missing","MOD11 (scaled)")),
|
|
| 314 | 334 |
rectangles=list(border=NA,col=c(NA,"tan","darkgreen")), |
| 315 | 335 |
text=list(c("C5 MOD35 Processing Path","Desert","Land")),
|
| 316 |
rectangles=list(border=NA,col=c(NA,IGBP$col[sort(unique(transd$value[transd$variable=="MCD12Q1"]+1))])),
|
|
| 336 |
rectangles=list(border=NA,col=c(NA,IGBP$col[sort(unique(transd$value[transd$variable=="MCD12Q1"]+1))])), |
|
| 317 | 337 |
text=list(c("MCD12Q1 IGBP Land Cover",IGBP$class[sort(unique(transd$value[transd$variable=="MCD12Q1"]+1))])))))),
|
| 318 | 338 |
strip = strip.custom(par.strip.text=list(cex=.75))) |
| 319 | 339 |
print(p4) |
Also available in: Unified diff
Added script to test functionality of HEG tool which reveals that v2.12 segfaults when multiple bands are selected for processing. Also updated MOD35_ExtractProcessPath to process one-band-at-a-time instead of in batches.