Revision 4f171339
Added by Adam M. Wilson about 11 years ago
| climate/procedures/MOD35_Explore.r | ||
|---|---|---|
| 1 | 1 |
## explore the MOD35 data downloaded and gridded by the DAAC |
| 2 | 2 |
setwd("~/acrobates/adamw/projects/interp/data/modis/mod35")
|
| 3 |
|
|
| 3 |
setwd("/Users/adamw/GoogleDrive/Work/presentations/20130426_NASAMeeting/data")
|
|
| 4 | 4 |
library(raster) |
| 5 | 5 |
library(rasterVis) |
| 6 | 6 |
library(rgdal) |
| ... | ... | |
| 28 | 28 |
|
| 29 | 29 |
## get % cloudy |
| 30 | 30 |
v5=stack(brick("../mod06/summary/MOD06_h11v08_ymoncld01.nc",varname="CLD01"))
|
| 31 |
v5=stack(brick("summary/MOD06_h11v08_ymoncld01.nc",varname="CLD01"))
|
|
| 31 | 32 |
projection(v5)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs" |
| 32 | 33 |
v6=stack(brick("summary/MOD35_h11v08.nc",varname="PCloud"))
|
| 33 | 34 |
projection(v6)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs" |
| ... | ... | |
| 39 | 40 |
v5m=mean(v5) |
| 40 | 41 |
|
| 41 | 42 |
|
| 43 |
## worldclim |
|
| 44 |
wc=brick("summary/worldclim_h11v08.tif")
|
|
| 45 |
names(wc)=paste("WC",1:12,sep="_")
|
|
| 46 |
projection(wc)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs" |
|
| 47 |
wcsp=projectRaster(wc,crs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0") |
|
| 48 |
wcsp=readAll(wcsp) |
|
| 49 |
wcsp=setZ(wcsp,paste("2011-",1:12,"-15",sep=""))
|
|
| 50 |
|
|
| 42 | 51 |
## landcover |
| 43 | 52 |
lulc=raster("~/acrobatesroot/jetzlab/Data/environ/global/landcover/MODIS/MCD12Q1_IGBP_2005_v51.tif")
|
| 53 |
lulc=raster("~/acrobatesroot/jetzlab/Data/environ/global/landcover/MODIS/MCD12Q1_IGBP_2005_v51.tif")
|
|
| 54 |
|
|
| 44 | 55 |
projection(lulc)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs" |
| 45 | 56 |
lulc=crop(lulc,v6) |
| 46 | 57 |
|
| ... | ... | |
| 63 | 74 |
|
| 64 | 75 |
### load WORLDCLIM elevation |
| 65 | 76 |
dem=raster("../../tiles/h11v08/dem_h11v08.tif",format="GTiff")
|
| 77 |
|
|
| 66 | 78 |
projection(dem)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs" |
| 67 | 79 |
|
| 80 |
dem=raster("summary/dem_h11v08.tif",format="GTiff")
|
|
| 68 | 81 |
dif=v6-v5 |
| 69 | 82 |
names(dif)=month.name |
| 70 | 83 |
|
| ... | ... | |
| 83 | 96 |
### compare MOD43 and MOD17 products |
| 84 | 97 |
|
| 85 | 98 |
## MOD17 |
| 86 |
mod17=raster("../MOD17/Npp_1km_C5.1_mean_00_to_06.tif",format="GTiff")
|
|
| 87 |
mod17=crop(projectRaster(mod17,v6,method="bilinear"),v6) |
|
| 88 |
NAvalue(mod17)=32767 |
|
| 99 |
#mod17=raster("../MOD17/Npp_1km_C5.1_mean_00_to_06.tif",format="GTiff")
|
|
| 100 |
#mod17=crop(projectRaster(mod17,v6,method="bilinear"),v6)
|
|
| 101 |
#NAvalue(mod17)=32767
|
|
| 89 | 102 |
|
| 90 | 103 |
mod17qc=raster("../MOD17/Npp_QC_1km_C5.1_mean_00_to_06.tif",format="GTiff")
|
| 104 |
mod17qc=raster("summary/Npp_QC_1km_C5.1_mean_00_to_06.tif",format="GTiff")
|
|
| 105 |
|
|
| 91 | 106 |
mod17qc=crop(projectRaster(mod17qc,v6,method="bilinear"),v6) |
| 92 | 107 |
mod17qc[mod17qc<0|mod17qc>100]=NA |
| 93 | 108 |
|
| 94 | 109 |
## MOD43 via earth engine |
| 95 |
mod43=raster("../mod43/3b21aa90cc657523ff31e9559f36fb12.EVI_MEAN.tif",format="GTiff")
|
|
| 96 |
mod43=crop(projectRaster(mod43,v6,method="bilinear"),v6) |
|
| 110 |
#mod43=raster("../mod43/3b21aa90cc657523ff31e9559f36fb12.EVI_MEAN.tif",format="GTiff")
|
|
| 111 |
#mod43=crop(projectRaster(mod43,v6,method="bilinear"),v6)
|
|
| 97 | 112 |
|
| 98 | 113 |
mod43qc=raster("../mod43/3b21aa90cc657523ff31e9559f36fb12.Percent_Cloudy.tif",format="GTiff")
|
| 114 |
mod43qc=raster("summary/3b21aa90cc657523ff31e9559f36fb12.Percent_Cloudy.tif",format="GTiff")
|
|
| 99 | 115 |
mod43qc=crop(projectRaster(mod43qc,v6,method="bilinear"),v6) |
| 100 | 116 |
mod43qc[mod43qc<0|mod43qc>100]=NA |
| 101 | 117 |
|
| ... | ... | |
| 103 | 119 |
modprod=stack(mod17qc,mod43qc) |
| 104 | 120 |
names(modprod)=c("MOD17","MOD43")
|
| 105 | 121 |
|
| 106 |
### |
|
| 122 |
### worldclim and interpolation |
|
| 123 |
m7=raster("summary/h11v08_pred_2_8.nc",varname="ppt")
|
|
| 124 |
m7=readAll(m7) |
|
| 125 |
m7=m7*30 |
|
| 126 |
pred=stack(subset(wc,2),subset(m7,2)) |
|
| 127 |
names(pred)=c("WorldClim","MOD35 Precipitation")
|
|
| 107 | 128 |
|
| 108 | 129 |
n=100 |
| 109 | 130 |
at=seq(0,100,len=n) |
| 110 | 131 |
cols=grey(seq(0,1,len=n)) |
| 111 | 132 |
cols=rainbow(n) |
| 112 | 133 |
bgyr=colorRampPalette(c("blue","green","yellow","red"))
|
| 134 |
gbg=colorRampPalette(c(grey(.2),"grey","blue")) |
|
| 135 |
|
|
| 113 | 136 |
cols=bgyr(n) |
| 114 | 137 |
|
| 115 | 138 |
#levelplot(lulcf,margin=F,layers="LULC") |
| 139 |
compare=stack(v5,v6) |
|
| 116 | 140 |
|
| 117 | 141 |
m=3 |
| 118 | 142 |
mcompare=stack(subset(v5,m),subset(v6,m)) |
| 143 |
names(compare)=paste(rep(c("C5","C6"),each=12),1:12,sep="_")
|
|
| 119 | 144 |
|
| 120 | 145 |
mdiff=subset(v5,m)-subset(v6,m) |
| 121 | 146 |
names(mcompare)=c("Collection_5","Collection_6")
|
| 122 | 147 |
names(mdiff)=c("Collection_5-Collection_6")
|
| 123 | 148 |
|
| 149 |
library(Cairo) |
|
| 150 |
(trellis.par.get()) |
|
| 151 |
trellis.par.set(par.ylab.text=list(cex=5)) |
|
| 124 | 152 |
|
| 125 |
CairoPDF("output/mod35compare.pdf",width=11,height=8.5)
|
|
| 153 |
pdf("output/mod35compare.pdf",width=11,height=8.5, fontsize=18)
|
|
| 126 | 154 |
#CairoPNG("output/mod35compare_%d.png",units="in", width=11,height=8.5,pointsize=4000,dpi=1200,antialias="subpixel")
|
| 127 | 155 |
|
| 128 |
levelplot(mcompare,col.regions=cols,at=at,margin=F,sub="Frequency of MOD35 Clouds in March") |
|
| 156 |
levelplot(mcompare,col.regions=cols,at=at,margin=F,sub="Cloud Frequency (%) from MOD35 in March") |
|
| 157 |
|
|
| 158 |
## show all months c5 vs c6 |
|
| 159 |
levelplot(compare,col.regions=cols,at=at,margin=F,sub="C5 and C6 MOD35 Comparison", |
|
| 160 |
layout=c(12,2),as.table=T) |
|
| 161 |
|
|
| 129 | 162 |
#levelplot(dif,col.regions=bgyr(20),margin=F) |
| 130 | 163 |
levelplot(mdiff,col.regions=bgyr(100),at=seq(mdiff@data@min,mdiff@data@max,len=100),margin=F) |
| 131 | 164 |
|
| 132 | 165 |
|
| 133 |
boxplot(as.matrix(subset(dif,subset=1))~forest,varwidth=T,notch=T);abline(h=0) |
|
| 166 |
#boxplot(as.matrix(subset(dif,subset=1))~forest,varwidth=T,notch=T);abline(h=0)
|
|
| 134 | 167 |
|
| 135 | 168 |
|
| 136 | 169 |
levelplot(modprod,main="Missing Data (%) in MOD17 (NPP) and MOD43 (BRDF Reflectance)", |
| ... | ... | |
| 140 | 173 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at, |
| 141 | 174 |
xlim=c(-7300000,-6670000),ylim=c(0,600000)) |
| 142 | 175 |
|
| 143 |
levelplot(v5m,main="Missing Data (%) in MOD17 (NPP) and MOD43 (BRDF Reflectance)", |
|
| 144 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at, |
|
| 145 |
xlim=c(-7200000,-6670000),ylim=c(0,400000),margin=F) |
|
| 146 |
|
|
| 147 |
|
|
| 148 |
levelplot(subset(v5v6compare,1:2),main="Proportion Cloudy Days (%) in Collection 5 and 6 MOD35", |
|
| 176 |
levelplot(subset(v5v6compare,1:2),main="Mean Annual Proportion Cloudy Days (%) in Collection 5 and 6 MOD35", |
|
| 149 | 177 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at, |
| 150 | 178 |
margin=F) |
| 151 | 179 |
|
| 152 |
levelplot(subset(v5v6compare,1:2),main="Proportion Cloudy Days (%) in Collection 5 and 6 MOD35", |
|
| 180 |
levelplot(subset(v5v6compare,1:2),main="Mean Annual Proportion Cloudy Days (%) in Collection 5 and 6 MOD35",
|
|
| 153 | 181 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at, |
| 154 | 182 |
xlim=c(-7200000,-6670000),ylim=c(0,400000),margin=F) |
| 155 | 183 |
|
| ... | ... | |
| 157 | 185 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at, |
| 158 | 186 |
xlim=c(-7500000,-7200000),ylim=c(700000,1000000),margin=F) |
| 159 | 187 |
|
| 188 |
levelplot(pred,main="Comparison of WorldClim with Cloud-Informed interpolation", |
|
| 189 |
sub="Tile H11v08 (Venezuela)",col.regions=grey(seq(0,1,len=n)),at=quantile(as.numeric(as.matrix(pred)),na.rm=T,seq(0,1,len=n))) |
|
| 190 |
|
|
| 191 |
levelplot(pred,main="Comparison of WorldClim with Cloud-Informed interpolation", |
|
| 192 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=quantile(as.numeric(as.matrix(pred)),na.rm=T,seq(0,1,len=n))) |
|
| 193 |
|
|
| 194 |
levelplot(pred,main="Comparison of WorldClim with Cloud-Informed interpolation", |
|
| 195 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=seq(0,360,len=n)) |
|
| 196 |
|
|
| 197 |
p.strip <- list(cex=1.5, lines=2, fontface='bold') |
|
| 198 |
ckey <- list(labels=list(cex=2), height=1) |
|
| 199 |
x.scale <- list(cex=1.5, alternating=1) |
|
| 200 |
y.scale <- list(cex=1.5, alternating=1) |
|
| 201 |
|
|
| 202 |
levelplot(pred,main="Comparison of WorldClim with Cloud-Informed interpolation", |
|
| 203 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=c(seq(-2,150,len=n*.75),seq(150.1,360,len=n*.25)), |
|
| 204 |
xlim=c(-7700000,-6700000),ylim=c(350000,1100000) , scales=list(x=x.scale,y=y.scale),par.strip.text=p.strip,colorkey=ckey, |
|
| 205 |
layout=c(2,1)) |
|
| 206 |
|
|
| 160 | 207 |
|
| 161 | 208 |
dev.off() |
| 162 | 209 |
|
| ... | ... | |
| 184 | 231 |
|
| 185 | 232 |
|
| 186 | 233 |
#### export KML timeseries |
| 234 |
wd="/Users/adamw/GoogleDrive/Work/presentations/20130426_NASAMeeting/data/" |
|
| 187 | 235 |
library(plotKML) |
| 188 | 236 |
tile="h11v08" |
| 189 | 237 |
file=paste("summary/MOD35_",tile,".nc",sep="")
|
| 190 |
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=""))
|
|
| 238 |
system(paste("gdalwarp -overwrite -multi -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=""))
|
|
| 191 | 239 |
|
| 192 |
v6sp=brick(paste("MOD35_",tile,".tif",sep=""))
|
|
| 240 |
v6sp=brick(paste("summary/MOD35_",tile,".tif",sep=""))
|
|
| 193 | 241 |
v6sp=readAll(v6sp) |
| 194 | 242 |
|
| 195 |
## wasn't working with line below, perhaps Z should just be text? not date? |
|
| 196 |
v6sp=setZ(v6sp,as.Date(paste("2011-",1:12,"-15",sep="")))
|
|
| 243 |
v6sp=setZ(v6sp,paste("2011-",1:12,"-15",sep=""))
|
|
| 197 | 244 |
names(v6sp)=month.name |
| 198 | 245 |
|
| 199 |
kml_open("output/mod35.kml")
|
|
| 246 |
|
|
| 247 |
setwd(paste(wd,"/kml/mod35",sep="")) |
|
| 248 |
kmlf="mod35.kml" |
|
| 249 |
kml_open(file.name=kmlf) |
|
| 200 | 250 |
|
| 201 | 251 |
|
| 202 | 252 |
kml_layer.RasterBrick(v6sp, |
| 203 | 253 |
plot.legend = TRUE, dtime = "", tz = "GMT", |
| 204 |
z.lim = c(0,100),colour_scale = get("colour_scale_numeric", envir = plotKML.opts))
|
|
| 254 |
z.lim = c(0,100),colour_scale = gbg(20),
|
|
| 205 | 255 |
# home_url = get("home_url", envir = plotKML.opts),
|
| 206 | 256 |
# metadata = NULL, html.table = NULL, |
| 207 |
# altitudeMode = "clampToGround", balloon = FALSE,
|
|
| 257 |
altitudeMode = "clampToGround", balloon = FALSE,folder.name=kmlfl,legend.file="mod35_legend.png"
|
|
| 208 | 258 |
) |
| 209 | 259 |
|
| 210 |
logo = "http://static.tumblr.com/t0afs9f/KWTm94tpm/yale_logo.png" |
|
| 211 |
kml_screen(image.file = logo, position = "UL", sname = "YALE logo",size=c(.1,.1)) |
|
| 212 |
kml_close("mod35.kml")
|
|
| 213 |
kml_compress("mod35.kml",files=c(paste(month.name,".png",sep=""),"obj_legend.png"),zip="/usr/bin/zip")
|
|
| 260 |
#logo = "yale_logo.png" |
|
| 261 |
#kml_screen(image.file = logo, position = "UL", sname = "YALE logo",size=c(.1,.1)) |
|
| 262 |
kml_close(file.name=kmlf) |
|
| 263 |
#kml_compress(file.name=kmlf,files=c(paste(month.name,".png",sep=""),"obj_legend.png"),zip="/usr/bin/zip") |
|
| 264 |
|
|
| 265 |
### worldclim |
|
| 266 |
setwd(paste(wd,"/kml/wc",sep="")) |
|
| 267 |
kmlf="worldclim.kml" |
|
| 268 |
kml_open(file.name=kmlf) |
|
| 269 |
|
|
| 270 |
kml_layer.RasterBrick(wcsp, |
|
| 271 |
plot.legend = TRUE, dtime = "", tz = "GMT", |
|
| 272 |
z.lim = c(0,650),colour_scale = gbg(100), |
|
| 273 |
# home_url = get("home_url", envir = plotKML.opts),
|
|
| 274 |
# metadata = NULL, html.table = NULL, |
|
| 275 |
altitudeMode = "clampToGround", balloon = FALSE |
|
| 276 |
) |
|
| 277 |
|
|
| 278 |
kml_close(file.name=kmlf) |
|
| 279 |
#kml_compress(file.name=kmlf,files=c(paste("WC_",1:12,".png",sep=""),"obj_legend.png"),zip="/usr/bin/zip")
|
|
| 280 |
|
|
| 281 |
setwd(wd) |
|
Also available in: Unified diff
Old code