Project

General

Profile

« Previous | Next » 

Revision b3f5074e

Added by Adam Wilson over 11 years ago

Added initial code to make KML of cloud climatologies

View differences:

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