Revision b3f5074e
Added by Adam Wilson about 11 years ago
| climate/procedures/MOD35_Explore.r | ||
|---|---|---|
| 1 | 1 |
## explore the MOD35 data downloaded and gridded by the DAAC |
| 2 |
setwd("~/acrobates/projects/interp/data/modis/mod35")
|
|
| 2 |
setwd("~/acrobates/adamw/projects/interp/data/modis/mod35")
|
|
| 3 | 3 |
|
| 4 | 4 |
library(raster) |
| 5 | 5 |
library(rasterVis) |
| ... | ... | |
| 29 | 29 |
## get % cloudy |
| 30 | 30 |
v5=stack(brick("../mod06/summary/MOD06_h11v08_ymoncld01.nc",varname="CLD01"))
|
| 31 | 31 |
projection(v5)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs" |
| 32 |
v6=stack(brick("MOD35_h11v08.nc",varname="CLD01"))
|
|
| 32 |
v6=stack(brick("summary/MOD35_h11v08.nc",varname="PCloud"))
|
|
| 33 | 33 |
projection(v6)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs" |
| 34 |
v6=setZ(v6,as.Date(paste("2011-",1:12,"-15",sep="")))
|
|
| 35 |
names(v6)=month.name |
|
| 34 | 36 |
|
| 35 | 37 |
## generate means |
| 36 | 38 |
v6m=mean(v6) |
| ... | ... | |
| 38 | 40 |
|
| 39 | 41 |
|
| 40 | 42 |
## landcover |
| 41 |
lulc=raster("~/acrobatesroot/Data/environ/global/landcover/MODIS/MCD12Q1_IGBP_2005_v51.tif")
|
|
| 43 |
lulc=raster("~/acrobatesroot/jetzlab/Data/environ/global/landcover/MODIS/MCD12Q1_IGBP_2005_v51.tif")
|
|
| 42 | 44 |
projection(lulc)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs" |
| 43 | 45 |
lulc=crop(lulc,v6) |
| 44 | 46 |
|
| ... | ... | |
| 59 | 61 |
levels(lulcf)=list(data.frame(ID=0:16,LULC=lulc_levels,LULC2=lulc_levels2)) |
| 60 | 62 |
|
| 61 | 63 |
|
| 62 |
|
|
| 63 | 64 |
### load WORLDCLIM elevation |
| 64 | 65 |
dem=raster("../../tiles/h11v08/dem_h11v08.tif",format="GTiff")
|
| 65 | 66 |
projection(dem)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs" |
| ... | ... | |
| 68 | 69 |
names(dif)=month.name |
| 69 | 70 |
|
| 70 | 71 |
difm=v6m-v5m |
| 72 |
v5v6compare=stack(v5m,v6m,difm) |
|
| 73 |
names(v5v6compare)=c("Collection 5","Collection 6","Difference (C6-C5)")
|
|
| 71 | 74 |
|
| 72 | 75 |
tile=extent(v6) |
| 73 | 76 |
|
| ... | ... | |
| 76 | 79 |
forest=lulcm>=1&lulcm<=5 |
| 77 | 80 |
|
| 78 | 81 |
|
| 79 |
boxplot(cld) |
|
| 80 |
splom(cld) |
|
| 81 |
|
|
| 82 |
|
|
| 83 | 82 |
##################################### |
| 84 | 83 |
### compare MOD43 and MOD17 products |
| 85 | 84 |
|
| ... | ... | |
| 104 | 103 |
modprod=stack(mod17qc,mod43qc) |
| 105 | 104 |
names(modprod)=c("MOD17","MOD43")
|
| 106 | 105 |
|
| 106 |
### |
|
| 107 |
|
|
| 107 | 108 |
n=100 |
| 108 | 109 |
at=seq(0,100,len=n) |
| 109 | 110 |
cols=grey(seq(0,1,len=n)) |
| ... | ... | |
| 121 | 122 |
names(mdiff)=c("Collection_5-Collection_6")
|
| 122 | 123 |
|
| 123 | 124 |
|
| 124 |
pdf("output/mod35compare.pdf",width=11,height=8.5)
|
|
| 125 |
CairoPDF("output/mod35compare.pdf",width=11,height=8.5)
|
|
| 126 |
#CairoPNG("output/mod35compare_%d.png",units="in", width=11,height=8.5,pointsize=4000,dpi=1200,antialias="subpixel")
|
|
| 125 | 127 |
|
| 126 | 128 |
levelplot(mcompare,col.regions=cols,at=at,margin=F,sub="Frequency of MOD35 Clouds in March") |
| 127 |
levelplot(dif,col.regions=bgyr(20),margin=F) |
|
| 128 |
levelplot(mdiff,col.regions=bgyr(20),margin=F)
|
|
| 129 |
#levelplot(dif,col.regions=bgyr(20),margin=F)
|
|
| 130 |
levelplot(mdiff,col.regions=bgyr(100),at=seq(mdiff@data@min,mdiff@data@max,len=100),margin=F)
|
|
| 129 | 131 |
|
| 130 | 132 |
|
| 131 | 133 |
boxplot(as.matrix(subset(dif,subset=1))~forest,varwidth=T,notch=T);abline(h=0) |
| 132 | 134 |
|
| 133 |
dev.off() |
|
| 134 |
|
|
| 135 | 135 |
|
| 136 | 136 |
levelplot(modprod,main="Missing Data (%) in MOD17 (NPP) and MOD43 (BRDF Reflectance)", |
| 137 | 137 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at) |
| ... | ... | |
| 144 | 144 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at, |
| 145 | 145 |
xlim=c(-7200000,-6670000),ylim=c(0,400000),margin=F) |
| 146 | 146 |
|
| 147 |
levelplot(stack(v5m,v6m),main="Missing Data (%) in MOD17 (NPP) and MOD43 (BRDF Reflectance)", |
|
| 147 |
|
|
| 148 |
levelplot(subset(v5v6compare,1:2),main="Proportion Cloudy Days (%) in Collection 5 and 6 MOD35", |
|
| 149 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at, |
|
| 150 |
margin=F) |
|
| 151 |
|
|
| 152 |
levelplot(subset(v5v6compare,1:2),main="Proportion Cloudy Days (%) in Collection 5 and 6 MOD35", |
|
| 148 | 153 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at, |
| 149 | 154 |
xlim=c(-7200000,-6670000),ylim=c(0,400000),margin=F) |
| 150 | 155 |
|
| 156 |
levelplot(subset(v5v6compare,1:2),main="Proportion Cloudy Days (%) in Collection 5 and 6 MOD35", |
|
| 157 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at, |
|
| 158 |
xlim=c(-7500000,-7200000),ylim=c(700000,1000000),margin=F) |
|
| 159 |
|
|
| 160 |
|
|
| 161 |
dev.off() |
|
| 162 |
|
|
| 151 | 163 |
### smoothing plots |
| 152 | 164 |
## explore smoothed version |
| 153 | 165 |
td=subset(v6,m) |
| ... | ... | |
| 168 | 180 |
|
| 169 | 181 |
62,59 |
| 170 | 182 |
0,3 |
| 183 |
|
|
| 184 |
|
|
| 185 |
|
|
| 186 |
#### export KML timeseries |
|
| 187 |
library(plotKML) |
|
| 188 |
tile="h11v08" |
|
| 189 |
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=""))
|
|
| 191 |
|
|
| 192 |
v6sp=brick(paste("MOD35_",tile,".tif",sep=""))
|
|
| 193 |
v6sp=readAll(v6sp) |
|
| 194 |
|
|
| 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="")))
|
|
| 197 |
names(v6sp)=month.name |
|
| 198 |
|
|
| 199 |
kml_open("output/mod35.kml")
|
|
| 200 |
|
|
| 201 |
|
|
| 202 |
kml_layer.RasterBrick(v6sp, |
|
| 203 |
plot.legend = TRUE, dtime = "", tz = "GMT", |
|
| 204 |
z.lim = c(0,100),colour_scale = get("colour_scale_numeric", envir = plotKML.opts))
|
|
| 205 |
# home_url = get("home_url", envir = plotKML.opts),
|
|
| 206 |
# metadata = NULL, html.table = NULL, |
|
| 207 |
# altitudeMode = "clampToGround", balloon = FALSE, |
|
| 208 |
) |
|
| 209 |
|
|
| 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")
|
|
Also available in: Unified diff
Added initial code to make KML of cloud climatologies