| 1 |
ddd9a810
|
Adam M. Wilson
|
## Figures associated with MOD35 Cloud Mask Exploration
|
| 2 |
|
|
|
| 3 |
|
|
setwd("~/acrobates/adamw/projects/MOD35C5")
|
| 4 |
|
|
|
| 5 |
|
|
library(raster);beginCluster(10)
|
| 6 |
|
|
library(rasterVis)
|
| 7 |
|
|
library(rgdal)
|
| 8 |
|
|
library(plotKML)
|
| 9 |
|
|
library(Cairo)
|
| 10 |
|
|
library(reshape)
|
| 11 |
|
|
|
| 12 |
|
|
## get % cloudy
|
| 13 |
|
|
mod09=raster("data/MOD09_2009.tif")
|
| 14 |
|
|
names(mod09)="MOD09_cloud"
|
| 15 |
|
|
|
| 16 |
|
|
mod35c5=raster("data/MOD35_2009.tif")
|
| 17 |
|
|
mod35c5=crop(mod35c5,mod09)
|
| 18 |
|
|
names(mod35c5)="MOD35C5_cloud"
|
| 19 |
|
|
|
| 20 |
|
|
mod35c6=raster("")
|
| 21 |
|
|
|
| 22 |
|
|
## landcover
|
| 23 |
|
|
if(!file.exists("data/MCD12Q1_IGBP_2005_v51_1km_wgs84.tif")){
|
| 24 |
|
|
system(paste("gdalwarp -r near -ot Byte -co \"COMPRESS=LZW\"",
|
| 25 |
|
|
" ~/acrobatesroot/jetzlab/Data/environ/global/landcover/MODIS/MCD12Q1_IGBP_2005_v51.tif ",
|
| 26 |
|
|
" -t_srs \"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs\" ",
|
| 27 |
|
|
tempdir(),"/MCD12Q1_IGBP_2005_v51_wgs84.tif -overwrite ",sep=""))
|
| 28 |
|
|
lulc=raster(paste(tempdir(),"/MCD12Q1_IGBP_2005_v51_wgs84.tif",sep=""))
|
| 29 |
|
|
## aggregate to 1km resolution
|
| 30 |
|
|
lulc2=aggregate(lulc,2,fun=function(x,na.rm=T) {
|
| 31 |
|
|
x=na.omit(x)
|
| 32 |
|
|
ux <- unique(x)
|
| 33 |
|
|
ux[which.max(tabulate(match(x, ux)))]
|
| 34 |
|
|
},file=paste(tempdir(),"/1km.tif",sep=""))
|
| 35 |
|
|
writeRaster(lulc2,"data/MCD12Q1_IGBP_2005_v51_1km_wgs84.tif",options=c("COMPRESS=LZW","ZLEVEL=9","PREDICTOR=2"),datatype="INT1U",overwrite=T)
|
| 36 |
|
|
}
|
| 37 |
|
|
lulc=raster("data/MCD12Q1_IGBP_2005_v51_1km_wgs84.tif")
|
| 38 |
|
|
# lulc=ratify(lulc)
|
| 39 |
|
|
data(worldgrids_pal) #load palette
|
| 40 |
|
|
IGBP=data.frame(ID=0:16,col=worldgrids_pal$IGBP[-c(18,19)],
|
| 41 |
|
|
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)
|
| 42 |
|
|
IGBP$class=rownames(IGBP);rownames(IGBP)=1:nrow(IGBP)
|
| 43 |
|
|
levels(lulc)=list(IGBP)
|
| 44 |
|
|
extent(lulc)=alignExtent(lulc,mod09)
|
| 45 |
|
|
names(lulc)="MCD12Q1"
|
| 46 |
|
|
|
| 47 |
|
|
## make land mask
|
| 48 |
|
|
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)
|
| 49 |
|
|
land=raster("data/land.tif")
|
| 50 |
|
|
|
| 51 |
|
|
#####################################
|
| 52 |
|
|
### compare MOD43 and MOD17 products
|
| 53 |
|
|
|
| 54 |
|
|
## MOD17
|
| 55 |
|
|
mod17=raster("data/MOD17A3_Science_NPP_mean_00_12.tif",format="GTiff")
|
| 56 |
|
|
NAvalue(mod17)=65535
|
| 57 |
|
|
#extent(mod17)=alignExtent(mod17,mod09)
|
| 58 |
|
|
mod17=crop(mod17,mod09)
|
| 59 |
|
|
names(mod17)="MOD17"
|
| 60 |
|
|
|
| 61 |
|
|
mod17qc=raster("data/MOD17A3_Science_NPP_Qc_mean_00_12.tif",format="GTiff")
|
| 62 |
|
|
NAvalue(mod17qc)=255
|
| 63 |
|
|
#extent(mod17qc)=alignExtent(mod17qc,mod09)
|
| 64 |
|
|
mod17qc=crop(mod17qc,mod09)
|
| 65 |
|
|
names(mod17qc)="MOD17qc"
|
| 66 |
|
|
|
| 67 |
|
|
## MOD11 via earth engine
|
| 68 |
|
|
mod11=raster("data/MOD11_2009.tif",format="GTiff")
|
| 69 |
|
|
names(mod11)="MOD11"
|
| 70 |
|
|
mod11qc=raster("data/MOD11_Pmiss_2009.tif",format="GTiff")
|
| 71 |
|
|
names(mod11qc)="MOD11qc"
|
| 72 |
|
|
|
| 73 |
|
|
## MOD43 via earth engine
|
| 74 |
|
|
mod43=raster("data/mod43_2009.tif",format="GTiff")
|
| 75 |
|
|
mod43qc=raster("data/mod43_2009.tif",format="GTiff")
|
| 76 |
|
|
|
| 77 |
|
|
|
| 78 |
|
|
### Create some summary objects for plotting
|
| 79 |
|
|
#difm=v6m-v5m
|
| 80 |
|
|
#v5v6compare=stack(v5m,v6m,difm)
|
| 81 |
|
|
#names(v5v6compare)=c("Collection 5","Collection 6","Difference (C6-C5)")
|
| 82 |
|
|
|
| 83 |
|
|
### Processing path
|
| 84 |
|
|
pp=raster("data/MOD35_ProcessPath.tif")
|
| 85 |
|
|
extent(pp)=alignExtent(pp,mod09)
|
| 86 |
|
|
pp=crop(pp,mod09)
|
| 87 |
|
|
|
| 88 |
|
|
## Summary plot of mod17 and mod43
|
| 89 |
|
|
modprod=stack(mod35c5,mod09,pp,lulc)#,mod43,mod43qc)
|
| 90 |
|
|
names(modprod)=c("MOD17","MOD17qc")#,"MOD43","MOD43qc")
|
| 91 |
|
|
|
| 92 |
|
|
|
| 93 |
|
|
## comparison of % cloudy days
|
| 94 |
|
|
dif=mod35c5-mod09
|
| 95 |
|
|
hist(dif,maxsamp=1000000)
|
| 96 |
|
|
|
| 97 |
|
|
## draw lulc-stratified random sample of mod35-mod09 differences
|
| 98 |
|
|
samp=sampleStratified(lulc, 1000, exp=10)
|
| 99 |
|
|
save(samp,file="LULC_StratifiedSample_10000.Rdata")
|
| 100 |
|
|
|
| 101 |
|
|
mean(dif[samp],na.rm=T)
|
| 102 |
|
|
|
| 103 |
|
|
Stats(dif,function(x) c(mean=mean(x),sd=sd(x)))
|
| 104 |
|
|
|
| 105 |
|
|
|
| 106 |
|
|
###
|
| 107 |
|
|
|
| 108 |
|
|
n=100
|
| 109 |
|
|
at=seq(0,100,len=n)
|
| 110 |
|
|
cols=grey(seq(0,1,len=n))
|
| 111 |
|
|
cols=rainbow(n)
|
| 112 |
|
|
bgyr=colorRampPalette(c("blue","green","yellow","red"))
|
| 113 |
|
|
cols=bgyr(n)
|
| 114 |
|
|
|
| 115 |
|
|
#levelplot(lulcf,margin=F,layers="LULC")
|
| 116 |
|
|
|
| 117 |
|
|
CairoPDF("output/mod35compare.pdf",width=11,height=8.5)
|
| 118 |
|
|
#CairoPNG("output/mod35compare_%d.png",units="in", width=11,height=8.5,pointsize=4000,dpi=1200,antialias="subpixel")
|
| 119 |
|
|
|
| 120 |
|
|
### Transects
|
| 121 |
|
|
r1=Lines(list(
|
| 122 |
|
|
Line(matrix(c(
|
| 123 |
|
|
-61.183,1.165,
|
| 124 |
|
|
-60.881,0.825
|
| 125 |
|
|
),ncol=2,byrow=T))),"Venezuela")
|
| 126 |
|
|
r2=Lines(list(
|
| 127 |
|
|
Line(matrix(c(
|
| 128 |
|
|
133.746,-31.834,
|
| 129 |
|
|
134.226,-32.143
|
| 130 |
|
|
),ncol=2,byrow=T))),"Australia")
|
| 131 |
|
|
r3=Lines(list(
|
| 132 |
|
|
Line(matrix(c(
|
| 133 |
|
|
73.943,27.419,
|
| 134 |
|
|
74.369,26.877
|
| 135 |
|
|
),ncol=2,byrow=T))),"India")
|
| 136 |
|
|
r4=Lines(list(
|
| 137 |
|
|
Line(matrix(c(
|
| 138 |
|
|
-5.164,42.270,
|
| 139 |
|
|
-4.948,42.162
|
| 140 |
|
|
),ncol=2,byrow=T))),"Spain")
|
| 141 |
|
|
|
| 142 |
|
|
r5=Lines(list(
|
| 143 |
|
|
Line(matrix(c(
|
| 144 |
|
|
24.170,-17.769,
|
| 145 |
|
|
24.616,-18.084
|
| 146 |
|
|
),ncol=2,byrow=T))),"Africa")
|
| 147 |
|
|
|
| 148 |
|
|
|
| 149 |
|
|
trans=SpatialLines(list(r1,r2,r3,r4,r5),CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs "))
|
| 150 |
|
|
|
| 151 |
|
|
transd=lapply(list(mod35c5,mod09,mod17,mod17qc,mod11qc,lulc,pp),function(l) {
|
| 152 |
|
|
td=extract(l,trans,along=T,cellnumbers=F)
|
| 153 |
|
|
names(td)=names(trans) # colnames(td)=c("value","transect")
|
| 154 |
|
|
cells=extract(l,trans,along=T,cellnumbers=T)
|
| 155 |
|
|
cells2=lapply(cells,function(x) xyFromCell(l,x[,1]))
|
| 156 |
|
|
dists=lapply(cells2,function(x) spDistsN1(x,x[1,],longlat=T))
|
| 157 |
|
|
td2=do.call(rbind.data.frame,lapply(1:length(td),function(i) cbind.data.frame(value=td[[i]],cells2[[i]],dist=dists[[i]],transect=names(td)[i])))
|
| 158 |
|
|
td2$prod=names(l)
|
| 159 |
|
|
td2$loc=rownames(td2)
|
| 160 |
|
|
td2=td2[order(td2$dist),]
|
| 161 |
|
|
print(paste("Finished ",names(l)))
|
| 162 |
|
|
return(td2)}
|
| 163 |
|
|
)
|
| 164 |
|
|
transdl=melt(transd,id.vars=c("prod","transect","loc","x","y","dist"))
|
| 165 |
|
|
transd$loc=as.numeric(transd$loc)
|
| 166 |
|
|
transdl$type=ifelse(grepl("MOD35|MOD09|qc",transdl$prod),"QC","Data")
|
| 167 |
|
|
|
| 168 |
|
|
nppid=transdl$prod=="MOD17"
|
| 169 |
|
|
|
| 170 |
|
|
xyplot(value~dist|transect,groups=prod,type=c("smooth","p"),
|
| 171 |
|
|
data=transdl,panel=function(...,subscripts=subscripts) {
|
| 172 |
|
|
td=transdl[subscripts,]
|
| 173 |
|
|
## mod09
|
| 174 |
|
|
imod09=td$prod=="MOD09_cloud"
|
| 175 |
|
|
panel.xyplot(td$dist[imod09],td$value[imod09],type=c("p","smooth"),span=0.2,subscripts=1:sum(imod09),col="red",pch=16,cex=.5)
|
| 176 |
|
|
## mod35C5
|
| 177 |
|
|
imod35=td$prod=="MOD35C5_cloud"
|
| 178 |
|
|
panel.xyplot(td$dist[imod35],td$value[imod35],type=c("p","smooth"),span=0.09,subscripts=1:sum(imod35),col="blue",pch=16,cex=.5)
|
| 179 |
|
|
## mod17
|
| 180 |
|
|
imod17=td$prod=="MOD17"
|
| 181 |
|
|
panel.xyplot(td$dist[imod17],100*td$value[imod17]/max(td$value[imod17]),
|
| 182 |
|
|
type=c("smooth"),span=0.09,subscripts=1:sum(imod17),col="darkgreen",lty="dashed",pch=1,cex=.5)
|
| 183 |
|
|
imod17qc=td$prod=="MOD17qc"
|
| 184 |
|
|
panel.xyplot(td$dist[imod17qc],td$value[imod17qc],type=c("p","smooth"),span=0.09,subscripts=1:sum(imod17qc),col="darkgreen",pch=16,cex=.5)
|
| 185 |
|
|
## mod11
|
| 186 |
|
|
# imod11=td$prod=="MOD11"
|
| 187 |
|
|
# panel.xyplot(td$dist[imod17],100*td$value[imod17]/max(td$value[imod17]),
|
| 188 |
|
|
# type=c("smooth"),span=0.09,subscripts=1:sum(imod17),col="darkgreen",lty="dashed",pch=1,cex=.5)
|
| 189 |
|
|
imod11qc=td$prod=="MOD11qc"
|
| 190 |
|
|
panel.xyplot(td$dist[imod11qc],td$value[imod11qc],type=c("p","smooth"),span=0.09,subscripts=1:sum(imod11qc),col="maroon",pch=16,cex=.5)
|
| 191 |
|
|
## means
|
| 192 |
|
|
means=td$prod%in%c("","MOD17qc","MOD09_cloud","MOD35_cloud")
|
| 193 |
|
|
## land
|
| 194 |
|
|
path=td[td$prod=="MOD35_ProcessPath",]
|
| 195 |
|
|
panel.segments(path$dist,0,c(path$dist[-1],max(path$dist)),0,col=IGBP$col[path$value],subscripts=1:nrow(path),lwd=15,type="l")
|
| 196 |
|
|
land=td[td$prod=="MCD12Q1",]
|
| 197 |
|
|
panel.segments(land$dist,-5,c(land$dist[-1],max(land$dist)),-5,col=IGBP$col[land$value],subscripts=1:nrow(land),lwd=15,type="l")
|
| 198 |
|
|
},subscripts=T,par.settings = list(grid.pars = list(lineend = "butt")),
|
| 199 |
|
|
scales=list(
|
| 200 |
|
|
x=list(alternating=1), #lim=c(0,50),
|
| 201 |
|
|
y=list(at=c(-5,0,seq(20,100,len=5)),
|
| 202 |
|
|
labels=c("IGBP","MOD35",seq(20,100,len=5)),
|
| 203 |
|
|
lim=c(-10,100))),
|
| 204 |
|
|
xlab="Distance Along Transect (km)",
|
| 205 |
|
|
key=list(space="right",lines=list(col=c("red","blue","darkgreen","maroon")),text=list(c("MOD09 % Cloudy","MOD35 % Cloudy","MOD17 % Missing","MOD11 % Missing"),lwd=1,col=c("red","blue","darkgreen","maroon"))))
|
| 206 |
|
|
|
| 207 |
|
|
### levelplot of regions
|
| 208 |
|
|
|
| 209 |
|
|
|
| 210 |
|
|
c(levelplot(mod35c5,margin=F),levelplot(mod09,margin=F),levelplot(mod11qc),levelplot(mod17qc),x.same = T, y.same = T)
|
| 211 |
|
|
|
| 212 |
|
|
levelplot(modprod)
|
| 213 |
|
|
|
| 214 |
|
|
|
| 215 |
|
|
|
| 216 |
|
|
### LANDCOVER
|
| 217 |
|
|
levelplot(lulcf,col.regions=levels(lulcf)[[1]]$col,
|
| 218 |
|
|
scales=list(cex=2),
|
| 219 |
|
|
colorkey=list(space="right",at=0:16,labels=list(at=seq(0.5,16.5,by=1),labels=levels(lulcf)[[1]]$class,cex=2)),margin=F)
|
| 220 |
|
|
|
| 221 |
|
|
|
| 222 |
|
|
levelplot(mcompare,col.regions=cols,at=at,margin=F,sub="Frequency of MOD35 Clouds in March")
|
| 223 |
|
|
#levelplot(dif,col.regions=bgyr(20),margin=F)
|
| 224 |
|
|
levelplot(mdiff,col.regions=bgyr(100),at=seq(mdiff@data@min,mdiff@data@max,len=100),margin=F)
|
| 225 |
|
|
|
| 226 |
|
|
|
| 227 |
|
|
boxplot(as.matrix(subset(dif,subset=1))~forest,varwidth=T,notch=T);abline(h=0)
|
| 228 |
|
|
|
| 229 |
|
|
|
| 230 |
|
|
levelplot(modprod,main="Missing Data (%) in MOD17 (NPP) and MOD43 (BRDF Reflectance)",
|
| 231 |
|
|
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at)
|
| 232 |
|
|
|
| 233 |
|
|
|
| 234 |
|
|
|
| 235 |
|
|
|
| 236 |
|
|
levelplot(modprod,main="Missing Data (%) in MOD17 (NPP) and MOD43 (BRDF Reflectance)",
|
| 237 |
|
|
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at,
|
| 238 |
|
|
xlim=c(-7300000,-6670000),ylim=c(0,600000))
|
| 239 |
|
|
|
| 240 |
|
|
levelplot(v5m,main="Missing Data (%) in MOD17 (NPP) and MOD43 (BRDF Reflectance)",
|
| 241 |
|
|
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at,
|
| 242 |
|
|
xlim=c(-7200000,-6670000),ylim=c(0,400000),margin=F)
|
| 243 |
|
|
|
| 244 |
|
|
|
| 245 |
|
|
levelplot(subset(v5v6compare,1:2),main="Proportion Cloudy Days (%) in Collection 5 and 6 MOD35",
|
| 246 |
|
|
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at,
|
| 247 |
|
|
margin=F)
|
| 248 |
|
|
|
| 249 |
|
|
levelplot(subset(v5v6compare,1:2),main="Proportion Cloudy Days (%) in Collection 5 and 6 MOD35",
|
| 250 |
|
|
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at,
|
| 251 |
|
|
xlim=c(-7200000,-6670000),ylim=c(0,400000),margin=F)
|
| 252 |
|
|
|
| 253 |
|
|
levelplot(subset(v5v6compare,1:2),main="Proportion Cloudy Days (%) in Collection 5 and 6 MOD35",
|
| 254 |
|
|
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at,
|
| 255 |
|
|
xlim=c(-7500000,-7200000),ylim=c(700000,1000000),margin=F)
|
| 256 |
|
|
|
| 257 |
|
|
|
| 258 |
|
|
dev.off()
|
| 259 |
|
|
|
| 260 |
|
|
### smoothing plots
|
| 261 |
|
|
## explore smoothed version
|
| 262 |
|
|
td=subset(v6,m)
|
| 263 |
|
|
## build weight matrix
|
| 264 |
|
|
s=3
|
| 265 |
|
|
w=matrix(1/(s*s),nrow=s,ncol=s)
|
| 266 |
|
|
#w[s-1,s-1]=4/12; w
|
| 267 |
|
|
td2=focal(td,w=w)
|
| 268 |
|
|
td3=stack(td,td2)
|
| 269 |
|
|
|
| 270 |
|
|
levelplot(td3,col.regions=cols,at=at,margin=F)
|
| 271 |
|
|
|
| 272 |
|
|
dev.off()
|
| 273 |
|
|
plot(stack(difm,lulc))
|
| 274 |
|
|
|
| 275 |
|
|
### ROI
|
| 276 |
|
|
tile_ll=projectExtent(v6, "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
|
| 277 |
|
|
|
| 278 |
|
|
62,59
|
| 279 |
|
|
0,3
|
| 280 |
|
|
|
| 281 |
|
|
|
| 282 |
|
|
|
| 283 |
|
|
#### export KML timeseries
|
| 284 |
|
|
library(plotKML)
|
| 285 |
|
|
tile="h11v08"
|
| 286 |
|
|
file=paste("summary/MOD35_",tile,".nc",sep="")
|
| 287 |
|
|
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=""))
|
| 288 |
|
|
|
| 289 |
|
|
v6sp=brick(paste("MOD35_",tile,".tif",sep=""))
|
| 290 |
|
|
v6sp=readAll(v6sp)
|
| 291 |
|
|
|
| 292 |
|
|
## wasn't working with line below, perhaps Z should just be text? not date?
|
| 293 |
|
|
v6sp=setZ(v6sp,as.Date(paste("2011-",1:12,"-15",sep="")))
|
| 294 |
|
|
names(v6sp)=month.name
|
| 295 |
|
|
|
| 296 |
|
|
kml_open("output/mod35.kml")
|
| 297 |
|
|
|
| 298 |
|
|
|
| 299 |
|
|
kml_layer.RasterBrick(v6sp,
|
| 300 |
|
|
plot.legend = TRUE, dtime = "", tz = "GMT",
|
| 301 |
|
|
z.lim = c(0,100),colour_scale = get("colour_scale_numeric", envir = plotKML.opts))
|
| 302 |
|
|
# home_url = get("home_url", envir = plotKML.opts),
|
| 303 |
|
|
# metadata = NULL, html.table = NULL,
|
| 304 |
|
|
# altitudeMode = "clampToGround", balloon = FALSE,
|
| 305 |
|
|
)
|
| 306 |
|
|
|
| 307 |
|
|
logo = "http://static.tumblr.com/t0afs9f/KWTm94tpm/yale_logo.png"
|
| 308 |
|
|
kml_screen(image.file = logo, position = "UL", sname = "YALE logo",size=c(.1,.1))
|
| 309 |
|
|
kml_close("mod35.kml")
|
| 310 |
|
|
kml_compress("mod35.kml",files=c(paste(month.name,".png",sep=""),"obj_legend.png"),zip="/usr/bin/zip")
|
