Revision 2bca6aaa
Added by Adam Wilson about 11 years ago
| climate/procedures/NDP-026D.R | ||
|---|---|---|
| 8 | 8 |
library(rasterVis) |
| 9 | 9 |
library(rgdal) |
| 10 | 10 |
## register parallel processing |
| 11 |
registerDoMC(20)
|
|
| 11 |
registerDoMC(10)
|
|
| 12 | 12 |
|
| 13 | 13 |
|
| 14 | 14 |
## available here http://cdiac.ornl.gov/epubs/ndp/ndp026d/ndp026d.html |
| ... | ... | |
| 110 | 110 |
projection(cldys)=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
|
| 111 | 111 |
|
| 112 | 112 |
#### Evaluate MOD35 Cloud data |
| 113 |
mod35=brick("../modis/mod35/MOD35_h11v08.nc",varname="CLD01")
|
|
| 114 |
mod35sd=brick("../modis/mod35/MOD35_h11v08.nc",varname="CLD_sd")
|
|
| 115 |
projection(mod35)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs" |
|
| 113 |
mod35c6=brick("~/acrobates/adamw/projects/MOD35C5/data/MOD35C6_2009_new.tif")
|
|
| 114 |
#projection(mod35c6)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs" |
|
| 116 | 115 |
|
| 117 | 116 |
|
| 118 | 117 |
### use data from google earth engine |
| 119 |
mod35=raster("../modis/mod09/global_2009/MOD35_2009.tif")
|
|
| 118 |
mod35c5=raster("../modis/mod09/global_2009/MOD35_2009.tif")
|
|
| 120 | 119 |
mod09=raster("../modis/mod09/global_2009/MOD09_2009.tif")
|
| 121 | 120 |
|
| 122 | 121 |
## LULC |
| 123 |
system(paste("gdalwarp -r near -co \"COMPRESS=LZW\" -tr ",paste(res(mod09),collapse=" ",sep=""),
|
|
| 124 |
"-tap -multi -t_srs \"", projection(mod09),"\" /mnt/data/jetzlab/Data/environ/global/landcover/MODIS/MCD12Q1_IGBP_2005_v51.tif ../modis/mod12/MCD12Q1_IGBP_2005_v51.tif")) |
|
| 122 |
#system(paste("gdalwarp -r near -co \"COMPRESS=LZW\" -tr ",paste(res(mod09),collapse=" ",sep=""),
|
|
| 123 |
# "-tap -multi -t_srs \"", projection(mod09),"\" /mnt/data/jetzlab/Data/environ/global/landcover/MODIS/MCD12Q1_IGBP_2005_v51.tif ../modis/mod12/MCD12Q1_IGBP_2005_v51.tif"))
|
|
| 125 | 124 |
lulc=raster("../modis/mod12/MCD12Q1_IGBP_2005_v51.tif")
|
| 126 | 125 |
lulc=ratify(lulc) |
| 127 | 126 |
require(plotKML); data(worldgrids_pal) #load IGBP palette |
| ... | ... | |
| 135 | 134 |
colr=colorRampPalette(c("black","green","red"))
|
| 136 | 135 |
cols=colr(n) |
| 137 | 136 |
|
| 138 |
dif=mod35-mod09 |
|
| 139 |
bwplot(dif~as.factor(lulc)) |
|
| 137 |
#dif=mod35-mod09
|
|
| 138 |
#bwplot(dif~as.factor(lulc))
|
|
| 140 | 139 |
|
| 141 |
levelplot(mod35,col.regions=cols,at=at,margins=F,maxpixels=1e6)#,xlim=c(-100,-50),ylim=c(0,10)) |
|
| 142 |
levelplot(lulc,att="class",col.regions=levels(lulc)[[1]]$col,margin=F,maxpixels=1e6) |
|
| 140 |
#levelplot(mod35,col.regions=cols,at=at,margins=F,maxpixels=1e6)#,xlim=c(-100,-50),ylim=c(0,10))
|
|
| 141 |
#levelplot(lulc,att="class",col.regions=levels(lulc)[[1]]$col,margin=F,maxpixels=1e6)
|
|
| 143 | 142 |
|
| 144 | 143 |
#cldys=spTransform(cldys,CRS(projection(mod35))) |
| 145 | 144 |
|
| 146 |
mod35v=foreach(m=unique(cldm$month),.combine="rbind") %do% {
|
|
| 147 |
dr=subset(mod35,subset=m);projection(dr)=projection(mod35) |
|
| 148 |
dr2=subset(mod35sd,subset=m);projection(dr2)=projection(mod35) |
|
| 149 |
ds=cldms[cldms$month==m,] |
|
| 150 |
ds$mod35=unlist(extract(dr,ds,buffer=10,fun=mean,na.rm=T)) |
|
| 145 |
#mod35v=foreach(m=unique(cldm$month),.combine="rbind") %do% {
|
|
| 146 |
# dr=subset(mod35,subset=m);projection(dr)=projection(mod35)
|
|
| 147 |
# dr2=subset(mod35sd,subset=m);projection(dr2)=projection(mod35)
|
|
| 148 |
# ds=cldms[cldms$month==m,]
|
|
| 149 |
# ds$mod35=unlist(extract(dr,ds,buffer=10,fun=mean,na.rm=T))
|
|
| 151 | 150 |
# ds$mod35sd=extract(dr2,ds,buffer=10) |
| 152 |
print(m) |
|
| 153 |
return(ds@data[!is.na(ds$mod35),])} |
|
| 151 |
# print(m)
|
|
| 152 |
# return(ds@data[!is.na(ds$mod35),])}
|
|
| 154 | 153 |
|
| 155 | 154 |
y=2009 |
| 156 | 155 |
d=cldys[cldys$year==y,] |
| 157 | 156 |
|
| 158 |
d$mod35_10=unlist(extract(mod35,d,buffer=10000,fun=mean,na.rm=T)) |
|
| 157 |
d$mod35c6_10=unlist(extract(mod35c6,d,buffer=10000,fun=mean,na.rm=T)) |
|
| 158 |
d$mod35c5_10=unlist(extract(mod35c5,d,buffer=10000,fun=mean,na.rm=T)) |
|
| 159 | 159 |
d$mod09_10=unlist(extract(mod09,d,buffer=10000,fun=mean,na.rm=T)) |
| 160 |
d$dif=d$mod35_10-d$mod09_10 |
|
| 161 |
d$dif2=d$mod35_10-d$cld |
|
| 160 |
#d$dif=d$mod35_10-d$mod09_10
|
|
| 161 |
#d$dif2=d$mod35_10-d$cld
|
|
| 162 | 162 |
|
| 163 | 163 |
d$lulc=unlist(extract(lulc,d)) |
| 164 | 164 |
d$lulc_10=unlist(extract(lulc,d,buffer=10000,fun=mode,na.rm=T)) |
| ... | ... | |
| 175 | 175 |
|
| 176 | 176 |
### exploratory plots |
| 177 | 177 |
xyplot(cld~mod09_10,groups=lulc,data=d@data,pch=16,cex=.5)+layer(panel.abline(0,1,col="red")) |
| 178 |
xyplot(cld~mod09_10+mod35_10|as.factor(lulc),data=d@data,type=c("p","r"),pch=16,cex=.25,auto.key=T)+layer(panel.abline(0,1,col="green"))
|
|
| 178 |
xyplot(cld~mod09_10+mod35c5_10|as.factor(lulc),data=d@data,type=c("p","r"),pch=16,cex=.25,auto.key=T)+layer(panel.abline(0,1,col="green"))
|
|
| 179 | 179 |
xyplot(cld~mod35_10|as.factor(lulc),data=d@data,pch=16,cex=.5)+layer(panel.abline(0,1,col="red")) |
| 180 | 180 |
xyplot(mod35_10~mod09_10|as.factor(lulc),data=d@data,pch=16,cex=.5)+layer(panel.abline(0,1,col="red")) |
| 181 | 181 |
|
Also available in: Unified diff
Updated cloud climatology evaluation