Revision 84d4d760
Added by Adam Wilson about 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 `find /home/adamw/acrobates/adamw/projects/interp/data/modis/mod35/summary/ -name '*h[1]*_mean.nc'` ")
|
|
| 26 |
system("gdalbuildvrt data/MOD35C6.vrt `find /home/adamw/acrobates/adamw/projects/interp/data/modis/mod35/summary/ -name '*h[1]*_mean.nc'` ")
|
|
| 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-9]*_mean.nc'` ")
|
|
| 26 |
# system("gdalbuildvrt data/MOD35C6.vrt `find /home/adamw/acrobates/adamw/projects/interp/data/modis/mod35/summary/ -name '*h[1-9]*_mean.nc'` ")
|
|
| 27 | 27 |
|
| 28 |
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 4 -b 1 data/MOD35C6_CFday_pmiss.vrt `find /home/adamw/acrobates/adamw/projects/interp/data/modis/mod35/summary/ -name '*h[1]*.nc'` ")
|
|
| 29 |
system("gdalwarp data/MOD35C6_CFday_pmiss.vrt data/MOD35C6_CFday_pmiss.tif -r bilinear")
|
|
| 28 |
# 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 4 -b 1 data/MOD35C6_CFday_pmiss.vrt `find /home/adamw/acrobates/adamw/projects/interp/data/modis/mod35/summary/ -name '*h[1]*.nc'` ")
|
|
| 29 |
# system("gdalwarp data/MOD35C6_CFday_pmiss.vrt data/MOD35C6_CFday_pmiss.tif -r bilinear")
|
|
| 30 | 30 |
|
| 31 | 31 |
system("align.sh data/MOD35C6.vrt data/MOD09_2009.tif data/MOD35C6_2009.tif")
|
| 32 |
system("align.sh data/MOD35C6_CFday_pmiss.vrt data/MOD09_2009.tif data/MOD35C6_CFday_pmiss.tif")
|
|
| 32 |
# system("align.sh data/MOD35C6_CFday_pmiss.vrt data/MOD09_2009.tif data/MOD35C6_CFday_pmiss.tif")
|
|
| 33 | 33 |
} |
| 34 | 34 |
mod35c6=raster("data/MOD35C6_2009.tif")
|
| 35 | 35 |
names(mod35c6)="C6MOD35CF" |
| ... | ... | |
| 212 | 212 |
#dif_c6_09=mod35c6-mod09 |
| 213 | 213 |
#dif_c5_c6=mod35c5-mod35c6 |
| 214 | 214 |
|
| 215 |
## exploring various ways to compare cloud products along landcover or processing path edges |
|
| 216 |
t1=trd1[[1]] |
|
| 217 |
dif_p=calc(trd1[[1]], function(x) (x[1]-x[3])/(1-x[1])) |
|
| 218 |
edge=calc(edge(subset(t1,"MCD12Q1"),classes=T,type="inner"),function(x) ifelse(x==1,1,NA)) |
|
| 219 |
edgeb=buffer(edge,width=5000) |
|
| 220 |
edgeb=calc(edgeb,function(x) ifelse(is.na(x),0,1)) |
|
| 221 |
names(edge)="edge" |
|
| 222 |
names(edgeb)="edgeb" |
|
| 223 |
|
|
| 224 |
pedge=calc(edge(subset(t1,"MOD35pp"),classes=T,type="inner"),function(x) ifelse(x==1,1,NA)) |
|
| 225 |
pedgeb=buffer(pedge,width=3000) |
|
| 226 |
pedgeb=calc(pedgeb,function(x) ifelse(is.na(x),0,1)) |
|
| 227 |
names(pedgeb)="pedgeb" |
|
| 228 |
|
|
| 229 |
td1=as.data.frame(stack(t1,edgeb,pedgeb)) |
|
| 230 |
td1l=melt(td1,id.vars=c("pedgeb","edgeb","MOD35pp","MCD12Q1"),measure.vars=c("C5MOD35CF","C6MOD35CF","C5MOD09CF"))
|
|
| 231 |
td1l=td1l[td1l$pedgeb==1|td1l$edgeb==1,] |
|
| 232 |
|
|
| 233 |
cast(MOD35pp~MCD12Q1~variable,fun.aggregate="mean",data=td1l) |
|
| 234 |
|
|
| 235 |
## Moving window |
|
| 236 |
tiles=expand.grid(xmin=seq(-180,170,by=10),ymin=seq(-60,80,by=10)) |
|
| 237 |
tiles$xmax=tiles$xmin+10;tiles$ymax=tiles$ymin+10 |
|
| 238 |
|
|
| 239 |
|
|
| 240 |
registerDoMC(10) |
|
| 215 |
################################################################################## |
|
| 216 |
## Identify problematic areas with hard edges in cloud frequency |
|
| 241 | 217 |
############################ |
| 242 |
writeLines(c(""), "log.txt")
|
|
| 218 |
library(multicore) |
|
| 219 |
|
|
| 220 |
## set up processing chunks |
|
| 243 | 221 |
nrw=nrow(mod35c5) |
| 244 | 222 |
nby=10 |
| 245 | 223 |
nrwg=seq(1,nrw,by=nby) |
| 246 |
writeLines(paste("Processing ",length(nrwg)," groups"))
|
|
| 224 |
writeLines(paste("Processing ",length(nrwg)," groups and",nrw,"lines"))
|
|
| 247 | 225 |
|
| 226 |
## Parallel loop to conduct moving window analysis and quantify pixels with significant shifts across pp or lulc boundaries |
|
| 248 | 227 |
output=mclapply(nrwg,function(ti){
|
| 249 | 228 |
## Extract focal areas |
| 250 | 229 |
ngb=5 |
| 251 | 230 |
vals=getValuesFocal(mod35c5,ngb=ngb,row=ti,nrows=nby) |
| 231 |
vals_mod09=getValuesFocal(mod09,ngb=ngb,row=ti,nrows=nby) |
|
| 252 | 232 |
pp_ind=getValuesFocal(pp,ngb=ngb,row=ti,nrows=nby) |
| 253 | 233 |
lulc_ind=getValuesFocal(lulc,ngb=ngb,row=ti,nrows=nby) |
| 254 | 234 |
## processing path |
| ... | ... | |
| 258 | 238 |
tind2=tind1[!is.na(tind1)&!is.na(tval1)] #which classes exist without NAs? |
| 259 | 239 |
if(length(unique(tind2))<2) return(255) #if only one class, return 255 |
| 260 | 240 |
if(sort(table(tind2),dec=T)[2]<5) return(254) # if too few observations of class 2, return 254 |
| 261 |
round(kruskal.test(tval1,tind1)$p.value*100) # if it works, return p.value*100
|
|
| 241 |
return(round(kruskal.test(tval1,tind1)$p.value*100)) # if it works, return p.value*100
|
|
| 262 | 242 |
})),nrow=nby,ncol=ncol(mod35c5),byrow=T)) # turn it back into a raster |
| 263 | 243 |
## udpate raster and write it |
| 264 | 244 |
extent(pp_bias)=extent(mod35c5[ti:(ti+nby-1),1:ncol(mod35c5),drop=F]) |
| 265 | 245 |
projection(pp_bias)=projection(mod35c5) |
| 246 |
NAvalue(pp_bias) <- 255 |
|
| 266 | 247 |
writeRaster(pp_bias,file=paste("data/tiles/pp_bias_",ti,".tif",sep=""),
|
| 267 | 248 |
format="GTiff",dataType="INT1U",overwrite=T,NAflag=255) #,options=c("COMPRESS=LZW","ZLEVEL=9")
|
| 268 | 249 |
## landcover |
| ... | ... | |
| 272 | 253 |
tind2=tind1[!is.na(tind1)&!is.na(tval1)] #which classes exist without NAs? |
| 273 | 254 |
if(length(unique(tind2))<2) return(255) #if only one class, return 255 |
| 274 | 255 |
if(sort(table(tind2),dec=T)[2]<5) return(254) # if too few observations of class 2, return 254 |
| 275 |
round(kruskal.test(tval1,tind1)$p.value*100) # if it works, return p.value*100
|
|
| 256 |
return(round(kruskal.test(tval1,tind1)$p.value*100)) # if it works, get p.value*100
|
|
| 276 | 257 |
})),nrow=nby,ncol=ncol(mod35c5),byrow=T)) # turn it back into a raster |
| 277 | 258 |
## udpate raster and write it |
| 278 | 259 |
extent(lulc_bias)=extent(mod35c5[ti:(ti+nby-1),1:ncol(mod35c5),drop=F]) |
| 279 | 260 |
projection(lulc_bias)=projection(mod35c5) |
| 261 |
NAvalue(lulc_bias) <- 255 |
|
| 280 | 262 |
writeRaster(lulc_bias,file=paste("data/tiles/lulc_bias_",ti,".tif",sep=""),
|
| 281 |
format="GTiff",dataType="INT1U",overwrite=T,NAflag=255) #,options=c("COMPRESS=LZW","ZLEVEL=9")
|
|
| 263 |
format="GTiff",dataType="INT1U",overwrite=T,NAflag=255)#,options=c("COMPRESS=LZW","ZLEVEL=9")
|
|
| 264 |
|
|
| 265 |
## MOD09 |
|
| 266 |
mod09_lulc_bias=raster(matrix(do.call(rbind,lapply(1:nrow(vals_mod09),function(i) {
|
|
| 267 |
tind1=lulc_ind[i,] #vector of indices |
|
| 268 |
tval1=vals_mod09[i,] # vector of values |
|
| 269 |
tind2=tind1[!is.na(tind1)&!is.na(tval1)] #which classes exist without NAs? |
|
| 270 |
if(length(unique(tind2))<2) return(255) #if only one class, return 255 |
|
| 271 |
if(sort(table(tind2),dec=T)[2]<5) return(254) # if too few observations of class 2, return 254 |
|
| 272 |
return(round(kruskal.test(tval1,tind1)$p.value*100)) # if it works, get p.value*100 |
|
| 273 |
})),nrow=nby,ncol=ncol(mod35c5),byrow=T)) # turn it back into a raster |
|
| 274 |
## udpate raster and write it |
|
| 275 |
extent(mod09_lulc_bias)=extent(mod09[ti:(ti+nby-1),1:ncol(mod09),drop=F]) |
|
| 276 |
projection(mod09_lulc_bias)=projection(mod09) |
|
| 277 |
NAvalue(mod09_lulc_bias) <- 255 |
|
| 278 |
writeRaster(mod09_lulc_bias,file=paste("data/tiles/mod09_lulc_bias_",ti,".tif",sep=""),
|
|
| 279 |
format="GTiff",dataType="INT1U",overwrite=T,NAflag=255)#,options=c("COMPRESS=LZW","ZLEVEL=9")
|
|
| 280 |
|
|
| 282 | 281 |
return(ti) |
| 283 |
},mc.cores=10) |
|
| 284 |
|
|
| 285 |
|
|
| 286 |
|
|
| 287 |
## original solution |
|
| 288 |
# pp_bias=raster(matrix(do.call(rbind,lapply(1:nrow(vals),function(i) {
|
|
| 289 |
# if(length(unique(na.omit(pp_ind[i,])))<2) return(255) |
|
| 290 |
# if(sort(table(pp_ind[i,]),dec=T)[2]<5) return(254) |
|
| 291 |
# kruskal.test(vals[i,],pp_ind[i,])$p.value*100 |
|
| 292 |
# })),nrow=nrow(t_mod35c5),ncol=ncol(t_mod35c5),byrow=T)) |
|
| 293 |
# extent(pp_bias)=extent(t_mod35c5) |
|
| 294 |
# writeRaster(pp_bias,file=paste("data/tiles/pp_bias_",ti,".tif",sep=""),
|
|
| 295 |
# format="GTiff",dataType="INT1U",options=c("COMPRESS=LZW","ZLEVEL=9"),overwrite=T)
|
|
| 296 |
## Run kruskal test for processing lulc bias |
|
| 297 |
# lulc_bias=raster(matrix(do.call(rbind,mclapply(1:nrow(vals),function(i) {
|
|
| 298 |
# if(length(unique(na.omit(lulc_ind[i,])))<2) return(NA) |
|
| 299 |
# if(sort(table(lulc_ind[i,]),dec=T)[2]<5) return(255) |
|
| 300 |
# kruskal.test(vals[i,],lulc_ind[i,])$p.value*100 |
|
| 301 |
# })),nrow=nrow(t_mod35c5),ncol=ncol(t_mod35c5),byrow=T)) |
|
| 302 |
# extent(lulc_bias)=extent(t_mod35c5) |
|
| 303 |
# writeRaster(lulc_bias,dataType="INT1U",file=paste("data/tiles/lulc_bias_",ti,".tif",sep=""),
|
|
| 304 |
# format="GTiff",options=c("COMPRESS=LZW","ZLEVEL=9"),overwrite=T)
|
|
| 282 |
},mc.cores=15) |
|
| 305 | 283 |
|
| 306 | 284 |
|
| 307 | 285 |
### check raster temporary files |
| ... | ... | |
| 309 | 287 |
#removeTmpFiles(h=1) |
| 310 | 288 |
|
| 311 | 289 |
## merge all the files back again |
| 312 |
system("gdalbuildvrt data/lulc_bias.vrt `find data/tiles -name 'lulc_bias*tif'` ")
|
|
| 313 |
system("gdalwarp -co 'COMPRESS=LZW' -co 'ZLEVEL=9' data/lulc_bias.vrt data/lulc_bias.tif -r near")
|
|
| 290 |
system("gdalbuildvrt -srcnodata 255 -vrtnodata 255 data/lulc_bias.vrt `find data/tiles -name 'lulc_bias*tif'` ")
|
|
| 291 |
system("gdalwarp -srcnodata 255 -dstnodata 255 -multi -r bilinear -co 'COMPRESS=LZW' -co 'ZLEVEL=9' data/lulc_bias.vrt data/lulc_bias.tif -r near")
|
|
| 292 |
# system("align.sh data/lulc_bias.vrt data/MOD09_2009.tif data/lulc_bias.tif")
|
|
| 314 | 293 |
|
| 315 | 294 |
system("gdalbuildvrt data/pp_bias.vrt `find data/tiles -name 'pp_bias*tif'` ")
|
| 316 |
system("gdalwarp -co 'COMPRESS=LZW' -co 'ZLEVEL=9' data/pp_bias.vrt data/pp_bias.tif -r near")
|
|
| 295 |
system("gdalwarp -srcnodata 255 -dstnodata 255 -multi -r bilinear -co 'COMPRESS=LZW' -co 'ZLEVEL=9' data/pp_bias.vrt data/pp_bias.tif -r near")
|
|
| 296 |
system("align.sh -srcnodata 255 -dstnodata 255 -multi -r bilinear data/pp_bias.vrt data/MOD09_2009.tif data/pp_bias_align.tif &")
|
|
| 317 | 297 |
|
| 298 |
system("gdalbuildvrt data/mod09_lulc_bias.vrt `find data/tiles -name 'mod09_lulc_bias*tif'` ")
|
|
| 299 |
system("gdalwarp -srcnodata 255 -dstnodata 255 -multi -r bilinear -co 'COMPRESS=LZW' -co 'ZLEVEL=9' data/mod09_lulc_bias.vrt data/mod09_lulc_bias.tif -r near")
|
|
| 318 | 300 |
|
| 319 |
#plot(stack(foc,x1)) |
|
| 320 |
pat=c(-0.02,seq(0,0.1,len=50),seq(0.1,1,len=50)) |
|
| 321 |
grayr2=colorRampPalette(c("red",grey(c(.75,.5,.25))))
|
|
| 322 |
levelplot(stack(pp_bias,lulc_bias),col.regions=c("cyan",grayr2(100)),at=pat,colorkey=list(at=pat,cuts=100),margin=F)
|
|
| 301 |
### read them back in |
|
| 302 |
pp_bias=raster("data/pp_bias.tif")
|
|
| 303 |
names(pp_bias)="Processing Path" |
|
| 304 |
lulc_bias=raster("data/lulc_bias.tif")
|
|
| 305 |
names(lulc_bias)="Land Use Land Cover" |
|
| 306 |
|
|
| 307 |
pat=c(0,0.05,1)#seq(0,0.-5,len=2) #,seq(0.05,.1,len=50)) |
|
| 308 |
grayr2=colorRampPalette(c("red","transparent"))#grey(c(.75,.5,.25))))
|
|
| 309 |
levelplot(stack(pp_bias,lulc_bias),col.regions=c(grayr2(2)),at=pat, |
|
| 310 |
colorkey=F,margin=F,maxpixels=1e6)+layer(sp.lines(coast,lwd=.5)) |
|
| 323 | 311 |
|
| 324 | 312 |
cor(td1$MOD17,td1$C6MOD35,use="complete",method="spearman") |
| 325 | 313 |
cor(td1$MOD17[td1$edgeb==1],td1$C5MOD35[td1$edgeb==1],use="complete",method="spearman") |
| 326 | 314 |
|
| 327 | 315 |
bwplot(value~MOD35pp|variable,data=td1l[td1l$pedgeb==1,],horizontal=F) |
| 328 | 316 |
|
| 329 |
|
|
| 317 |
crosstab(dif_c5_09,pp) |
|
| 330 | 318 |
### Correlations |
| 331 | 319 |
#trdw=cast(trd,trans+x+y~variable,value="value") |
| 332 | 320 |
#cor(trdw$MOD17,trdw$C5MOD35,use="complete",method="spearman") |
Also available in: Unified diff
Exploring evalution options for mod35c5