| 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")
|
Exploring evalution options for mod35c5