Revision 57d44dd9
Added by Adam Wilson about 11 years ago
| climate/procedures/MOD35C5_Evaluation.r | ||
|---|---|---|
| 31 | 31 |
system("align.sh data/MOD35C6.vrt data/MOD09_2009.tif data/MOD35C6_2009.tif")
|
| 32 | 32 |
system("align.sh data/MOD35C6_CFday_pmiss.vrt data/MOD09_2009.tif data/MOD35C6_CFday_pmiss.tif")
|
| 33 | 33 |
} |
| 34 |
mod35c6=raster("data/MOD35C6_2009_v1.tif")
|
|
| 34 |
mod35c6=raster("data/MOD35C6_2009.tif")
|
|
| 35 | 35 |
names(mod35c6)="C6MOD35CF" |
| 36 | 36 |
NAvalue(mod35c6)=255 |
| 37 | 37 |
|
| ... | ... | |
| 92 | 92 |
|
| 93 | 93 |
### Processing path |
| 94 | 94 |
if(!file.exists("data/MOD35pp.tif"))
|
| 95 |
system("align.sh data/MOD35_ProcessPath.tif data/MOD09_2009.tif data/MOD35pp.tif")
|
|
| 95 |
system("align.sh data/C5MOD35_ProcessPath.tif data/MOD09_2009.tif data/MOD35pp.tif")
|
|
| 96 | 96 |
pp=raster("data/MOD35pp.tif")
|
| 97 | 97 |
NAvalue(pp)=255 |
| 98 | 98 |
names(pp)="MOD35pp" |
| ... | ... | |
| 132 | 132 |
73.943,27.419, |
| 133 | 133 |
74.369,26.877 |
| 134 | 134 |
),ncol=2,byrow=T))),"India") |
| 135 |
#r4=Lines(list( |
|
| 136 |
# Line(matrix(c( |
|
| 137 |
# -5.164,42.270, |
|
| 138 |
# -4.948,42.162 |
|
| 139 |
# ),ncol=2,byrow=T))),"Spain") |
|
| 140 |
|
|
| 141 |
r5=Lines(list( |
|
| 135 |
r4=Lines(list( |
|
| 142 | 136 |
Line(matrix(c( |
| 143 | 137 |
33.195,12.512, |
| 144 | 138 |
33.802,12.894 |
| 145 | 139 |
),ncol=2,byrow=T))),"Sudan") |
| 146 | 140 |
|
| 147 |
#r6=Lines(list( |
|
| 148 |
# Line(matrix(c( |
|
| 149 |
# -63.353,-10.746, |
|
| 150 |
# -63.376,-9.310 |
|
| 151 |
# ),ncol=2,byrow=T))),"Brazil") |
|
| 152 | 141 |
|
| 153 | 142 |
|
| 154 |
trans=SpatialLines(list(r1,r2,r3,r5),CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs "))
|
|
| 143 |
trans=SpatialLines(list(r1,r2,r3,r4),CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs "))
|
|
| 155 | 144 |
### write out shapefiles of transects |
| 156 | 145 |
writeOGR(SpatialLinesDataFrame(trans,data=data.frame(ID=names(trans)),match.ID=F),"output",layer="transects",driver="ESRI Shapefile",overwrite=T) |
| 157 | 146 |
|
| ... | ... | |
| 224 | 213 |
#dif_c5_c6=mod35c5-mod35c6 |
| 225 | 214 |
|
| 226 | 215 |
## exploring various ways to compare cloud products along landcover or processing path edges |
| 227 |
#t1=trd1[[1]] |
|
| 228 |
#dif_p=calc(trd1[[1]], function(x) (x[1]-x[3])/(1-x[1])) |
|
| 229 |
#edge=calc(edge(subset(t1,"MCD12Q1"),classes=T,type="inner"),function(x) ifelse(x==1,1,NA)) |
|
| 230 |
#edgeb=buffer(edge,width=5000) |
|
| 231 |
#edgeb=calc(edgeb,function(x) ifelse(is.na(x),0,1)) |
|
| 232 |
#names(edge)="edge" |
|
| 233 |
#names(edgeb)="edgeb" |
|
| 234 |
#td1=as.data.frame(stack(t1,edge,edgeb)) |
|
| 235 |
#cor(td1$MOD17,td1$C6MOD35,use="complete",method="spearman") |
|
| 236 |
#cor(td1$MOD17[td1$edgeb==1],td1$C5MOD35[td1$edgeb==1],use="complete",method="spearman") |
|
| 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) |
|
| 241 |
############################ |
|
| 242 |
writeLines(c(""), "log.txt")
|
|
| 243 |
nrw=nrow(mod35c5) |
|
| 244 |
nby=10 |
|
| 245 |
nrwg=seq(1,nrw,by=nby) |
|
| 246 |
writeLines(paste("Processing ",length(nrwg)," groups"))
|
|
| 247 |
|
|
| 248 |
output=mclapply(nrwg,function(ti){
|
|
| 249 |
## Extract focal areas |
|
| 250 |
ngb=5 |
|
| 251 |
vals=getValuesFocal(mod35c5,ngb=ngb,row=ti,nrows=nby) |
|
| 252 |
pp_ind=getValuesFocal(pp,ngb=ngb,row=ti,nrows=nby) |
|
| 253 |
lulc_ind=getValuesFocal(lulc,ngb=ngb,row=ti,nrows=nby) |
|
| 254 |
## processing path |
|
| 255 |
pp_bias=raster(matrix(do.call(rbind,lapply(1:nrow(vals),function(i) {
|
|
| 256 |
tind1=pp_ind[i,] #vector of indices |
|
| 257 |
tval1=vals[i,] # vector of values |
|
| 258 |
tind2=tind1[!is.na(tind1)&!is.na(tval1)] #which classes exist without NAs? |
|
| 259 |
if(length(unique(tind2))<2) return(255) #if only one class, return 255 |
|
| 260 |
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 |
|
| 262 |
})),nrow=nby,ncol=ncol(mod35c5),byrow=T)) # turn it back into a raster |
|
| 263 |
## udpate raster and write it |
|
| 264 |
extent(pp_bias)=extent(mod35c5[ti:(ti+nby-1),1:ncol(mod35c5),drop=F]) |
|
| 265 |
projection(pp_bias)=projection(mod35c5) |
|
| 266 |
writeRaster(pp_bias,file=paste("data/tiles/pp_bias_",ti,".tif",sep=""),
|
|
| 267 |
format="GTiff",dataType="INT1U",overwrite=T,NAflag=255) #,options=c("COMPRESS=LZW","ZLEVEL=9")
|
|
| 268 |
## landcover |
|
| 269 |
lulc_bias=raster(matrix(do.call(rbind,lapply(1:nrow(vals),function(i) {
|
|
| 270 |
tind1=lulc_ind[i,] #vector of indices |
|
| 271 |
tval1=vals[i,] # vector of values |
|
| 272 |
tind2=tind1[!is.na(tind1)&!is.na(tval1)] #which classes exist without NAs? |
|
| 273 |
if(length(unique(tind2))<2) return(255) #if only one class, return 255 |
|
| 274 |
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 |
|
| 276 |
})),nrow=nby,ncol=ncol(mod35c5),byrow=T)) # turn it back into a raster |
|
| 277 |
## udpate raster and write it |
|
| 278 |
extent(lulc_bias)=extent(mod35c5[ti:(ti+nby-1),1:ncol(mod35c5),drop=F]) |
|
| 279 |
projection(lulc_bias)=projection(mod35c5) |
|
| 280 |
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")
|
|
| 282 |
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)
|
|
| 305 |
|
|
| 306 |
|
|
| 307 |
### check raster temporary files |
|
| 308 |
showTmpFiles() |
|
| 309 |
#removeTmpFiles(h=1) |
|
| 310 |
|
|
| 311 |
## 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")
|
|
| 314 |
|
|
| 315 |
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")
|
|
| 317 |
|
|
| 318 |
|
|
| 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)
|
|
| 323 |
|
|
| 324 |
cor(td1$MOD17,td1$C6MOD35,use="complete",method="spearman") |
|
| 325 |
cor(td1$MOD17[td1$edgeb==1],td1$C5MOD35[td1$edgeb==1],use="complete",method="spearman") |
|
| 326 |
|
|
| 327 |
bwplot(value~MOD35pp|variable,data=td1l[td1l$pedgeb==1,],horizontal=F) |
|
| 328 |
|
|
| 237 | 329 |
|
| 238 | 330 |
### Correlations |
| 239 | 331 |
#trdw=cast(trd,trans+x+y~variable,value="value") |
| ... | ... | |
| 252 | 344 |
## set up some graphing parameters |
| 253 | 345 |
at=seq(0,100,leng=100) |
| 254 | 346 |
bgyr=colorRampPalette(c("purple","blue","green","yellow","orange","red","red"))
|
| 255 |
bgrayr=colorRampPalette(c("purple","blue","grey","red","red"))
|
|
| 347 |
bgray=colorRampPalette(c("purple","blue","deepskyblue4"))
|
|
| 348 |
grayr=colorRampPalette(c("grey","red","darkred"))
|
|
| 349 |
bgrayr=colorRampPalette(c("darkblue","blue","grey","red","purple"))
|
|
| 350 |
|
|
| 256 | 351 |
cols=bgyr(100) |
| 257 | 352 |
|
| 353 |
strip=strip.custom(par.strip.text=list(cex=.7),bg="transparent") |
|
| 354 |
|
|
| 258 | 355 |
## global map |
| 259 | 356 |
library(maptools) |
| 260 | 357 |
coast=map2SpatialLines(map("world", interior=FALSE, plot=FALSE),proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
|
| 261 | 358 |
|
| 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)) |
|
| 359 |
g1=levelplot(stack(mod35c5,mod35c6,mod09),xlab=" ",scales=list(x=list(draw=F),y=list(alternating=1)), |
|
| 360 |
col.regions=cols,at=at,cuts=length(at),maxpixels=1e6, |
|
| 361 |
colorkey=list(at=at))+ |
|
| 362 |
# layer(sp.polygons(bbs,lwd=5,col="black"))+ |
|
| 363 |
layer(sp.lines(coast,lwd=.5))+ |
|
| 364 |
layer(sp.points(coordinates(bbs),col="black",cex=2,pch=13,lwd=2)) |
|
| 365 |
|
|
| 366 |
### Plot of differences between MOD09 adn MOD35 masks |
|
| 367 |
#system("gdalinfo -stats /home/adamw/acrobates/adamw/projects/MOD35C5/data/dif_c5_09.tif")
|
|
| 368 |
## get quantiles for color bar of differences |
|
| 369 |
#qs=unique(quantile(as.vector(as.matrix(dif_c5_09)),seq(0,1,len=100),na.rm=T)) |
|
| 370 |
#c(bgray(sum(qs<0)),grayr(sum(qs>=0)+1)) |
|
| 371 |
qs=seq(-80,80,len=100) |
|
| 372 |
g2=levelplot(dif_c5_09,col.regions=bgrayr(100),cuts=100,at=qs,margin=F,ylab=" ",colorkey=list("right",at=qs),maxpixels=1e6)+
|
|
| 373 |
layer(sp.points(coordinates(bbs),col="black",cex=2,pch=13,lwd=2))+ |
|
| 374 |
#layer(sp.polygons(bbs,lwd=2))+ |
|
| 375 |
layer(sp.lines(coast,lwd=.5)) |
|
| 263 | 376 |
|
| 264 |
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))
|
|
| 265 | 377 |
g2$strip=strip.custom(var.name="Difference (C5MOD35-C5MOD09)",style=1,strip.names=T,strip.levels=F) #update strip text |
| 266 | 378 |
#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)) |
| 267 | 379 |
|
| 268 | 380 |
### regional plots |
| 269 | 381 |
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"),
|
| 270 | 382 |
at=at,col.regions=cols,maxpixels=7e6, |
| 271 |
ylab="Latitude",xlab="Longitude"),strip = strip.custom(par.strip.text=list(cex=.7)))+layer(sp.lines(trans,lwd=2))
|
|
| 383 |
ylab="Latitude",xlab="Longitude"),strip.left=strip,strip = strip)+layer(sp.lines(trans,lwd=2))
|
|
| 272 | 384 |
|
| 273 | 385 |
p2=useOuterStrips( |
| 274 | 386 |
levelplot(value~x*y|variable+trans,data=trd[trd$variable%in%c("MCD12Q1"),],
|
| ... | ... | |
| 278 | 390 |
right=list(fun=draw.key(list(columns=1,#title="MCD12Q1 \n IGBP Land \n Cover", |
| 279 | 391 |
rectangles=list(col=IGBP$col,size=1), |
| 280 | 392 |
text=list(as.character(IGBP$ID),at=IGBP$ID-.5))))), |
| 281 |
ylab="",xlab=" "),strip = strip.custom(par.strip.text=list(cex=.7)),strip.left=F)+layer(sp.lines(trans,lwd=2))
|
|
| 393 |
ylab="",xlab=" "),strip = strip,strip.left=F)+layer(sp.lines(trans,lwd=2)) |
|
| 282 | 394 |
p3=useOuterStrips( |
| 283 | 395 |
levelplot(value~x*y|variable+trans,data=trd[trd$variable%in%c("MOD35pp"),],
|
| 284 | 396 |
asp=1,scales=list(draw=F,rot=0,relation="free"),colorkey=F, |
| ... | ... | |
| 287 | 399 |
right=list(fun=draw.key(list(columns=1,#title="MOD35 \n Processing \n Path", |
| 288 | 400 |
rectangles=list(col=c("blue","cyan","tan","darkgreen"),size=1),
|
| 289 | 401 |
text=list(c("Water","Coast","Desert","Land")))))),
|
| 290 |
ylab="",xlab=" "),strip = strip.custom(par.strip.text=list(cex=.7)),strip.left=F)+layer(sp.lines(trans,lwd=2))
|
|
| 402 |
ylab="",xlab=" "),strip = strip,strip.left=F)+layer(sp.lines(trans,lwd=2)) |
|
| 291 | 403 |
|
| 292 | 404 |
## transects |
| 293 | 405 |
p4=xyplot(value~dist|transect,groups=variable,type=c("smooth","p"),
|
| ... | ... | |
| 323 | 435 |
},subscripts=T,par.settings = list(grid.pars = list(lineend = "butt")), |
| 324 | 436 |
scales=list( |
| 325 | 437 |
x=list(alternating=1,relation="free"),#, lim=c(0,70)), |
| 326 |
y=list(at=c(-18,-10,seq(0,100,len=5)),
|
|
| 438 |
y=list(at=c(-20,-10,seq(0,100,len=5)),
|
|
| 327 | 439 |
labels=c("MCD12Q1 IGBP","MOD35 path",seq(0,100,len=5)),
|
| 328 | 440 |
lim=c(-25,100)), |
| 329 | 441 |
alternating=F), |
| ... | ... | |
| 339 | 451 |
text=list(c("C5 MOD35 Processing Path","Desert","Land")),
|
| 340 | 452 |
rectangles=list(border=NA,col=c(NA,IGBP$col[sort(unique(transd$value[transd$variable=="MCD12Q1"]+1))])), |
| 341 | 453 |
text=list(c("MCD12Q1 IGBP Land Cover",IGBP$class[sort(unique(transd$value[transd$variable=="MCD12Q1"]+1))])))))),
|
| 342 |
strip = strip.custom(par.strip.text=list(cex=.75))) |
|
| 343 |
print(p4) |
|
| 344 |
|
|
| 454 |
strip = strip,strip.left=F) |
|
| 455 |
#print(p4) |
|
| 345 | 456 |
|
| 346 | 457 |
|
| 347 | 458 |
CairoPDF("output/mod35compare.pdf",width=11,height=7)
|
| ... | ... | |
| 352 | 463 |
#print(g3,position=c(0.31,0.06,.42,0.27),more=F) |
| 353 | 464 |
|
| 354 | 465 |
### MOD35 Desert Processing path |
| 355 |
levelplot(pp,asp=1,scales=list(draw=T,rot=0),maxpixels=1e6, |
|
| 356 |
at=c(-1:3),col.regions=c("blue","cyan","tan","darkgreen"),margin=F,
|
|
| 357 |
colorkey=list(space="bottom",title="MOD35 Processing Path",labels=list(labels=c("Water","Coast","Desert","Land"),at=0:4-.5)))+layer(sp.polygons(bbs,lwd=2))+layer(sp.lines(coast,lwd=.5))
|
|
| 466 |
levelplot(pp,asp=1,scales=list(draw=T,rot=0),maxpixels=1e6,cuts=3, |
|
| 467 |
at=(0:3)+.5,col.regions=c("blue","cyan","tan","darkgreen"),margin=F,
|
|
| 468 |
colorkey=list(space="top",title="MOD35 Processing Path",labels=list(labels=c("Water","Coast","Desert","Land"),at=0:3),height=.25))+
|
|
| 469 |
layer(sp.points(coordinates(bbs),col="black",cex=2,pch=13,lwd=2))+ |
|
| 470 |
layer(sp.lines(coast,lwd=.5)) |
|
| 471 |
|
|
| 358 | 472 |
### levelplot of regions |
| 359 | 473 |
print(p1,position=c(0,0,.62,1),more=T) |
| 360 | 474 |
print(p2,position=c(0.6,0.21,0.78,0.79),more=T) |
| ... | ... | |
| 390 | 504 |
readAll(mod35c5) |
| 391 | 505 |
|
| 392 | 506 |
kml_layer.Raster(mod35c5, |
| 393 |
plot.legend = TRUE,raster_name="Collection 5 MOD35 Cloud Frequency", |
|
| 507 |
plot.legend = TRUE,raster_name="Collection 5 MOD35 2009 Cloud Frequency",
|
|
| 394 | 508 |
z.lim = c(0,100),colour_scale = get("colour_scale_numeric", envir = plotKML.opts),
|
| 395 | 509 |
# home_url = get("home_url", envir = plotKML.opts),
|
| 396 | 510 |
# metadata = NULL, html.table = NULL, |
Also available in: Unified diff
Added routine to process identify regions of concern in modis cloud data