Revision ddd9a810
Added by Adam Wilson over 11 years ago
| climate/procedures/EE_simple.py | ||
|---|---|---|
| 1 |
import ee |
|
| 2 |
from ee import mapclient |
|
| 3 |
|
|
| 4 |
import logging |
|
| 5 |
logging.basicConfig() |
|
| 6 |
|
|
| 7 |
MY_SERVICE_ACCOUNT = '511722844190@developer.gserviceaccount.com' # replace with your service account |
|
| 8 |
MY_PRIVATE_KEY_FILE = '/Users/adamw/EarthEngine-privatekey.p12' # replace with you private key file path |
|
| 9 |
|
|
| 10 |
ee.Initialize(ee.ServiceAccountCredentials(MY_SERVICE_ACCOUNT, MY_PRIVATE_KEY_FILE)) |
|
| 11 |
|
|
| 12 |
|
|
| 13 |
|
|
| 14 |
image = ee.Image('srtm90_v4')
|
|
| 15 |
map = ee.mapclient.MapClient() |
|
| 16 |
map.addOverlay(mapclient.MakeOverlay(image.getMapId({'min': 0, 'max': 3000})))
|
|
| climate/procedures/MOD35C5_Evaluation.r | ||
|---|---|---|
| 1 |
## 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")
|
|
| climate/procedures/MOD35C5_Evaluation_grass.r | ||
|---|---|---|
| 1 |
## 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 |
library(spgrass6) |
|
| 12 |
|
|
| 13 |
### set up grass database |
|
| 14 |
tf=paste(getwd(),"/grassdata",sep="") |
|
| 15 |
if(!file.exists(tf)) dir.create(tf) |
|
| 16 |
|
|
| 17 |
## set up GRASS |
|
| 18 |
gisBase="/usr/lib/grass64" |
|
| 19 |
initGRASS(gisBase=gisBase,override=T,gisDbase=tf,SG=as(raster("data/MOD17A3_Science_NPP_mean_00_12.tif"),SpatialGridDataFrame))
|
|
| 20 |
system(paste("g.proj -c proj4=\"",projection(glb),"\"",sep=""),ignore.stdout=T,ignore.stderr=T)
|
|
| 21 |
## read in NPP grid to serve as grid |
|
| 22 |
execGRASS("r.in.gdal",input=t@file@name,output="grid")
|
|
| 23 |
system("g.region rast=grid n=90 s=-90 save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)
|
|
| 24 |
|
|
| 25 |
|
|
| 26 |
## get % cloudy |
|
| 27 |
mod09=raster("data/MOD09_2009.tif")
|
|
| 28 |
names(mod09)="MOD09_cloud" |
|
| 29 |
|
|
| 30 |
mod35c5=raster("data/MOD35_2009.tif")
|
|
| 31 |
mod35c5=crop(mod35c5,mod09) |
|
| 32 |
names(mod35c5)="MOD35C5_cloud" |
|
| 33 |
|
|
| 34 |
mod35c6=raster("")
|
|
| 35 |
|
|
| 36 |
## landcover |
|
| 37 |
if(!file.exists("data/MCD12Q1_IGBP_2005_v51_1km_wgs84.tif")){
|
|
| 38 |
system(paste("gdalwarp -r near -ot Byte -co \"COMPRESS=LZW\"",
|
|
| 39 |
" ~/acrobatesroot/jetzlab/Data/environ/global/landcover/MODIS/MCD12Q1_IGBP_2005_v51.tif ", |
|
| 40 |
" -t_srs \"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs\" ", |
|
| 41 |
tempdir(),"/MCD12Q1_IGBP_2005_v51_wgs84.tif -overwrite ",sep="")) |
|
| 42 |
lulc=raster(paste(tempdir(),"/MCD12Q1_IGBP_2005_v51_wgs84.tif",sep="")) |
|
| 43 |
## aggregate to 1km resolution |
|
| 44 |
lulc2=aggregate(lulc,2,fun=function(x,na.rm=T) {
|
|
| 45 |
x=na.omit(x) |
|
| 46 |
ux <- unique(x) |
|
| 47 |
ux[which.max(tabulate(match(x, ux)))] |
|
| 48 |
},file=paste(tempdir(),"/1km.tif",sep="")) |
|
| 49 |
writeRaster(lulc2,"data/MCD12Q1_IGBP_2005_v51_1km_wgs84.tif",options=c("COMPRESS=LZW","ZLEVEL=9","PREDICTOR=2"),datatype="INT1U",overwrite=T)
|
|
| 50 |
} |
|
| 51 |
lulc=raster("data/MCD12Q1_IGBP_2005_v51_1km_wgs84.tif")
|
|
| 52 |
# lulc=ratify(lulc) |
|
| 53 |
data(worldgrids_pal) #load palette |
|
| 54 |
IGBP=data.frame(ID=0:16,col=worldgrids_pal$IGBP[-c(18,19)], |
|
| 55 |
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)
|
|
| 56 |
IGBP$class=rownames(IGBP);rownames(IGBP)=1:nrow(IGBP) |
|
| 57 |
levels(lulc)=list(IGBP) |
|
| 58 |
extent(lulc)=alignExtent(lulc,mod09) |
|
| 59 |
names(lulc)="MCD12Q1" |
|
| 60 |
|
|
| 61 |
## make land mask |
|
| 62 |
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)
|
|
| 63 |
land=raster("data/land.tif")
|
|
| 64 |
|
|
| 65 |
##################################### |
|
| 66 |
### compare MOD43 and MOD17 products |
|
| 67 |
|
|
| 68 |
## MOD17 |
|
| 69 |
mod17=raster("data/MOD17A3_Science_NPP_mean_00_12.tif",format="GTiff")
|
|
| 70 |
NAvalue(mod17)=65535 |
|
| 71 |
#extent(mod17)=alignExtent(mod17,mod09) |
|
| 72 |
mod17=crop(mod17,mod09) |
|
| 73 |
names(mod17)="MOD17" |
|
| 74 |
|
|
| 75 |
mod17qc=raster("data/MOD17A3_Science_NPP_Qc_mean_00_12.tif",format="GTiff")
|
|
| 76 |
NAvalue(mod17qc)=255 |
|
| 77 |
#extent(mod17qc)=alignExtent(mod17qc,mod09) |
|
| 78 |
mod17qc=crop(mod17qc,mod09) |
|
| 79 |
names(mod17qc)="MOD17qc" |
|
| 80 |
|
|
| 81 |
## MOD11 via earth engine |
|
| 82 |
mod11=raster("data/MOD11_2009.tif",format="GTiff")
|
|
| 83 |
names(mod11)="MOD11" |
|
| 84 |
mod11qc=raster("data/MOD11_Pmiss_2009.tif",format="GTiff")
|
|
| 85 |
names(mod11qc)="MOD11qc" |
|
| 86 |
|
|
| 87 |
## MOD43 via earth engine |
|
| 88 |
mod43=raster("data/mod43_2009.tif",format="GTiff")
|
|
| 89 |
mod43qc=raster("data/mod43_2009.tif",format="GTiff")
|
|
| 90 |
|
|
| 91 |
|
|
| 92 |
### Create some summary objects for plotting |
|
| 93 |
#difm=v6m-v5m |
|
| 94 |
#v5v6compare=stack(v5m,v6m,difm) |
|
| 95 |
#names(v5v6compare)=c("Collection 5","Collection 6","Difference (C6-C5)")
|
|
| 96 |
|
|
| 97 |
### Processing path |
|
| 98 |
pp=raster("data/MOD35_ProcessPath.tif")
|
|
| 99 |
extent(pp)=alignExtent(pp,mod09) |
|
| 100 |
pp=crop(pp,mod09) |
|
| 101 |
|
|
| 102 |
## Summary plot of mod17 and mod43 |
|
| 103 |
modprod=stack(mod35c5,mod09,pp,lulc)#,mod43,mod43qc) |
|
| 104 |
names(modprod)=c("MOD17","MOD17qc")#,"MOD43","MOD43qc")
|
|
| 105 |
|
|
| 106 |
|
|
| 107 |
## comparison of % cloudy days |
|
| 108 |
dif=mod35c5-mod09 |
|
| 109 |
hist(dif,maxsamp=1000000) |
|
| 110 |
|
|
| 111 |
## draw lulc-stratified random sample of mod35-mod09 differences |
|
| 112 |
samp=sampleStratified(lulc, 1000, exp=10) |
|
| 113 |
save(samp,file="LULC_StratifiedSample_10000.Rdata") |
|
| 114 |
|
|
| 115 |
mean(dif[samp],na.rm=T) |
|
| 116 |
|
|
| 117 |
Stats(dif,function(x) c(mean=mean(x),sd=sd(x))) |
|
| 118 |
|
|
| 119 |
|
|
| 120 |
### |
|
| 121 |
|
|
| 122 |
n=100 |
|
| 123 |
at=seq(0,100,len=n) |
|
| 124 |
cols=grey(seq(0,1,len=n)) |
|
| 125 |
cols=rainbow(n) |
|
| 126 |
bgyr=colorRampPalette(c("blue","green","yellow","red"))
|
|
| 127 |
cols=bgyr(n) |
|
| 128 |
|
|
| 129 |
#levelplot(lulcf,margin=F,layers="LULC") |
|
| 130 |
|
|
| 131 |
CairoPDF("output/mod35compare.pdf",width=11,height=8.5)
|
|
| 132 |
#CairoPNG("output/mod35compare_%d.png",units="in", width=11,height=8.5,pointsize=4000,dpi=1200,antialias="subpixel")
|
|
| 133 |
|
|
| 134 |
### Transects |
|
| 135 |
r1=Lines(list( |
|
| 136 |
Line(matrix(c( |
|
| 137 |
-61.183,1.165, |
|
| 138 |
-60.881,0.825 |
|
| 139 |
),ncol=2,byrow=T))),"Venezuela") |
|
| 140 |
r2=Lines(list( |
|
| 141 |
Line(matrix(c( |
|
| 142 |
133.746,-31.834, |
|
| 143 |
134.226,-32.143 |
|
| 144 |
),ncol=2,byrow=T))),"Australia") |
|
| 145 |
r3=Lines(list( |
|
| 146 |
Line(matrix(c( |
|
| 147 |
73.943,27.419, |
|
| 148 |
74.369,26.877 |
|
| 149 |
),ncol=2,byrow=T))),"India") |
|
| 150 |
r4=Lines(list( |
|
| 151 |
Line(matrix(c( |
|
| 152 |
-5.164,42.270, |
|
| 153 |
-4.948,42.162 |
|
| 154 |
),ncol=2,byrow=T))),"Spain") |
|
| 155 |
|
|
| 156 |
r5=Lines(list( |
|
| 157 |
Line(matrix(c( |
|
| 158 |
24.170,-17.769, |
|
| 159 |
24.616,-18.084 |
|
| 160 |
),ncol=2,byrow=T))),"Africa") |
|
| 161 |
|
|
| 162 |
|
|
| 163 |
trans=SpatialLines(list(r1,r2,r3,r4,r5),CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs "))
|
|
| 164 |
|
|
| 165 |
transd=lapply(list(mod35c5,mod09,mod17,mod17qc,mod11qc,lulc,pp),function(l) {
|
|
| 166 |
td=extract(l,trans,along=T,cellnumbers=F) |
|
| 167 |
names(td)=names(trans) # colnames(td)=c("value","transect")
|
|
| 168 |
cells=extract(l,trans,along=T,cellnumbers=T) |
|
| 169 |
cells2=lapply(cells,function(x) xyFromCell(l,x[,1])) |
|
| 170 |
dists=lapply(cells2,function(x) spDistsN1(x,x[1,],longlat=T)) |
|
| 171 |
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]))) |
|
| 172 |
td2$prod=names(l) |
|
| 173 |
td2$loc=rownames(td2) |
|
| 174 |
td2=td2[order(td2$dist),] |
|
| 175 |
print(paste("Finished ",names(l)))
|
|
| 176 |
return(td2)} |
|
| 177 |
) |
|
| 178 |
transdl=melt(transd,id.vars=c("prod","transect","loc","x","y","dist"))
|
|
| 179 |
transd$loc=as.numeric(transd$loc) |
|
| 180 |
transdl$type=ifelse(grepl("MOD35|MOD09|qc",transdl$prod),"QC","Data")
|
|
| 181 |
|
|
| 182 |
nppid=transdl$prod=="MOD17" |
|
| 183 |
|
|
| 184 |
xyplot(value~dist|transect,groups=prod,type=c("smooth","p"),
|
|
| 185 |
data=transdl,panel=function(...,subscripts=subscripts) {
|
|
| 186 |
td=transdl[subscripts,] |
|
| 187 |
## mod09 |
|
| 188 |
imod09=td$prod=="MOD09_cloud" |
|
| 189 |
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)
|
|
| 190 |
## mod35C5 |
|
| 191 |
imod35=td$prod=="MOD35C5_cloud" |
|
| 192 |
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)
|
|
| 193 |
## mod17 |
|
| 194 |
imod17=td$prod=="MOD17" |
|
| 195 |
panel.xyplot(td$dist[imod17],100*td$value[imod17]/max(td$value[imod17]), |
|
| 196 |
type=c("smooth"),span=0.09,subscripts=1:sum(imod17),col="darkgreen",lty="dashed",pch=1,cex=.5)
|
|
| 197 |
imod17qc=td$prod=="MOD17qc" |
|
| 198 |
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)
|
|
| 199 |
## mod11 |
|
| 200 |
# imod11=td$prod=="MOD11" |
|
| 201 |
# panel.xyplot(td$dist[imod17],100*td$value[imod17]/max(td$value[imod17]), |
|
| 202 |
# type=c("smooth"),span=0.09,subscripts=1:sum(imod17),col="darkgreen",lty="dashed",pch=1,cex=.5)
|
|
| 203 |
imod11qc=td$prod=="MOD11qc" |
|
| 204 |
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)
|
|
| 205 |
## means |
|
| 206 |
means=td$prod%in%c("","MOD17qc","MOD09_cloud","MOD35_cloud")
|
|
| 207 |
## land |
|
| 208 |
path=td[td$prod=="MOD35_ProcessPath",] |
|
| 209 |
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") |
|
| 210 |
land=td[td$prod=="MCD12Q1",] |
|
| 211 |
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") |
|
| 212 |
},subscripts=T,par.settings = list(grid.pars = list(lineend = "butt")), |
|
| 213 |
scales=list( |
|
| 214 |
x=list(alternating=1), #lim=c(0,50), |
|
| 215 |
y=list(at=c(-5,0,seq(20,100,len=5)), |
|
| 216 |
labels=c("IGBP","MOD35",seq(20,100,len=5)),
|
|
| 217 |
lim=c(-10,100))), |
|
| 218 |
xlab="Distance Along Transect (km)", |
|
| 219 |
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"))))
|
|
| 220 |
|
|
| 221 |
### levelplot of regions |
|
| 222 |
|
|
| 223 |
|
|
| 224 |
c(levelplot(mod35c5,margin=F),levelplot(mod09,margin=F),levelplot(mod11qc),levelplot(mod17qc),x.same = T, y.same = T) |
|
| 225 |
|
|
| 226 |
levelplot(modprod) |
|
| 227 |
|
|
| 228 |
|
|
| 229 |
|
|
| 230 |
### LANDCOVER |
|
| 231 |
levelplot(lulcf,col.regions=levels(lulcf)[[1]]$col, |
|
| 232 |
scales=list(cex=2), |
|
| 233 |
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) |
|
| 234 |
|
|
| 235 |
|
|
| 236 |
levelplot(mcompare,col.regions=cols,at=at,margin=F,sub="Frequency of MOD35 Clouds in March") |
|
| 237 |
#levelplot(dif,col.regions=bgyr(20),margin=F) |
|
| 238 |
levelplot(mdiff,col.regions=bgyr(100),at=seq(mdiff@data@min,mdiff@data@max,len=100),margin=F) |
|
| 239 |
|
|
| 240 |
|
|
| 241 |
boxplot(as.matrix(subset(dif,subset=1))~forest,varwidth=T,notch=T);abline(h=0) |
|
| 242 |
|
|
| 243 |
|
|
| 244 |
levelplot(modprod,main="Missing Data (%) in MOD17 (NPP) and MOD43 (BRDF Reflectance)", |
|
| 245 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at) |
|
| 246 |
|
|
| 247 |
|
|
| 248 |
|
|
| 249 |
|
|
| 250 |
levelplot(modprod,main="Missing Data (%) in MOD17 (NPP) and MOD43 (BRDF Reflectance)", |
|
| 251 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at, |
|
| 252 |
xlim=c(-7300000,-6670000),ylim=c(0,600000)) |
|
| 253 |
|
|
| 254 |
levelplot(v5m,main="Missing Data (%) in MOD17 (NPP) and MOD43 (BRDF Reflectance)", |
|
| 255 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at, |
|
| 256 |
xlim=c(-7200000,-6670000),ylim=c(0,400000),margin=F) |
|
| 257 |
|
|
| 258 |
|
|
| 259 |
levelplot(subset(v5v6compare,1:2),main="Proportion Cloudy Days (%) in Collection 5 and 6 MOD35", |
|
| 260 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at, |
|
| 261 |
margin=F) |
|
| 262 |
|
|
| 263 |
levelplot(subset(v5v6compare,1:2),main="Proportion Cloudy Days (%) in Collection 5 and 6 MOD35", |
|
| 264 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at, |
|
| 265 |
xlim=c(-7200000,-6670000),ylim=c(0,400000),margin=F) |
|
| 266 |
|
|
| 267 |
levelplot(subset(v5v6compare,1:2),main="Proportion Cloudy Days (%) in Collection 5 and 6 MOD35", |
|
| 268 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at, |
|
| 269 |
xlim=c(-7500000,-7200000),ylim=c(700000,1000000),margin=F) |
|
| 270 |
|
|
| 271 |
|
|
| 272 |
dev.off() |
|
| 273 |
|
|
| 274 |
### smoothing plots |
|
| 275 |
## explore smoothed version |
|
| 276 |
td=subset(v6,m) |
|
| 277 |
## build weight matrix |
|
| 278 |
s=3 |
|
| 279 |
w=matrix(1/(s*s),nrow=s,ncol=s) |
|
| 280 |
#w[s-1,s-1]=4/12; w |
|
| 281 |
td2=focal(td,w=w) |
|
| 282 |
td3=stack(td,td2) |
|
| 283 |
|
|
| 284 |
levelplot(td3,col.regions=cols,at=at,margin=F) |
|
| 285 |
|
|
| 286 |
dev.off() |
|
| 287 |
plot(stack(difm,lulc)) |
|
| 288 |
|
|
| 289 |
### ROI |
|
| 290 |
tile_ll=projectExtent(v6, "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") |
|
| 291 |
|
|
| 292 |
62,59 |
|
| 293 |
0,3 |
|
| 294 |
|
|
| 295 |
|
|
| 296 |
|
|
| 297 |
#### export KML timeseries |
|
| 298 |
library(plotKML) |
|
| 299 |
tile="h11v08" |
|
| 300 |
file=paste("summary/MOD35_",tile,".nc",sep="")
|
|
| 301 |
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=""))
|
|
| 302 |
|
|
| 303 |
v6sp=brick(paste("MOD35_",tile,".tif",sep=""))
|
|
| 304 |
v6sp=readAll(v6sp) |
|
| 305 |
|
|
| 306 |
## wasn't working with line below, perhaps Z should just be text? not date? |
|
| 307 |
v6sp=setZ(v6sp,as.Date(paste("2011-",1:12,"-15",sep="")))
|
|
| 308 |
names(v6sp)=month.name |
|
| 309 |
|
|
| 310 |
kml_open("output/mod35.kml")
|
|
| 311 |
|
|
| 312 |
|
|
| 313 |
kml_layer.RasterBrick(v6sp, |
|
| 314 |
plot.legend = TRUE, dtime = "", tz = "GMT", |
|
| 315 |
z.lim = c(0,100),colour_scale = get("colour_scale_numeric", envir = plotKML.opts))
|
|
| 316 |
# home_url = get("home_url", envir = plotKML.opts),
|
|
| 317 |
# metadata = NULL, html.table = NULL, |
|
| 318 |
# altitudeMode = "clampToGround", balloon = FALSE, |
|
| 319 |
) |
|
| 320 |
|
|
| 321 |
logo = "http://static.tumblr.com/t0afs9f/KWTm94tpm/yale_logo.png" |
|
| 322 |
kml_screen(image.file = logo, position = "UL", sname = "YALE logo",size=c(.1,.1)) |
|
| 323 |
kml_close("mod35.kml")
|
|
| 324 |
kml_compress("mod35.kml",files=c(paste(month.name,".png",sep=""),"obj_legend.png"),zip="/usr/bin/zip")
|
|
| climate/procedures/MOD35_Explore.r | ||
|---|---|---|
| 1 |
## explore the MOD35 data downloaded and gridded by the DAAC |
|
| 2 |
setwd("~/acrobates/adamw/projects/interp/data/modis/mod35")
|
|
| 3 |
|
|
| 4 |
library(raster) |
|
| 5 |
library(rasterVis) |
|
| 6 |
library(rgdal) |
|
| 7 |
library(plotKML) |
|
| 8 |
|
|
| 9 |
#f=list.files(pattern="*.hdf") |
|
| 10 |
|
|
| 11 |
#Sys.setenv(GEOL_AS_GCPS = "PARTIAL") |
|
| 12 |
## try swath-grid with gdal |
|
| 13 |
#GDALinfo(f[1]) |
|
| 14 |
#system(paste("gdalinfo",f[2]))
|
|
| 15 |
#GDALinfo("HDF4_EOS:EOS_SWATH:\"MOD35_L2.A2000100.1445.006.2012252024758.hdf\":mod35:Cloud_Mask")
|
|
| 16 |
#system("gdalinfo HDF4_EOS:EOS_SWATH:\"data/modis/mod35/MOD35_L2.A2000100.1445.006.2012252024758.hdf\":mod35:Cloud_Mask")
|
|
| 17 |
#system("gdalwarp -overwrite -geoloc -order 2 -r near -s_srs \"EPSG:4326\" HDF4_EOS:EOS_SWATH:\"MOD35_L2.A2000100.1445.006.2012252024758.hdf\":mod35:Cloud_Mask cloudmask.tif")
|
|
| 18 |
#system("gdalwarp -overwrite -r near -s_srs \"EPSG:4326\" HDF4_EOS:EOS_SWATH:\"MOD35_L2.A2000100.1445.006.2012252024758.hdf\":mod35:Cloud_Mask:1 cloudmask2.tif")
|
|
| 19 |
#r=raster(f[1]) |
|
| 20 |
#extent(r) |
|
| 21 |
#st=lapply(f[1:10],raster) |
|
| 22 |
#str=lapply(2:length(st),function(i) union(extent(st[[i-1]]),extent(st[[i]])))[[length(st)-1]] |
|
| 23 |
#str=union(extent(h11v08),str) |
|
| 24 |
#b1=brick(lapply(st,function(stt) {
|
|
| 25 |
# x=crop(alignExtent(stt,str),h11v08) |
|
| 26 |
# return(x) |
|
| 27 |
#})) |
|
| 28 |
#c=brick(f[1:10]) |
|
| 29 |
|
|
| 30 |
## get % cloudy |
|
| 31 |
v5=stack(brick("../mod06/summary/MOD06_h11v08_ymoncld01.nc",varname="CLD01"))
|
|
| 32 |
projection(v5)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs" |
|
| 33 |
v6=stack(brick("summary/MOD35_h11v08.nc",varname="PCloud"))
|
|
| 34 |
projection(v6)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs" |
|
| 35 |
v6=setZ(v6,as.Date(paste("2011-",1:12,"-15",sep="")))
|
|
| 36 |
names(v6)=month.name |
|
| 37 |
|
|
| 38 |
## generate means |
|
| 39 |
v6m=mean(v6) |
|
| 40 |
v5m=mean(v5) |
|
| 41 |
|
|
| 42 |
|
|
| 43 |
## landcover |
|
| 44 |
lulc=raster("~/acrobatesroot/jetzlab/Data/environ/global/landcover/MODIS/MCD12Q1_IGBP_2005_v51.tif")
|
|
| 45 |
projection(lulc)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs" |
|
| 46 |
lulc=crop(lulc,v6) |
|
| 47 |
|
|
| 48 |
Mode <- function(x,na.rm=T) { #get MODE
|
|
| 49 |
x=na.omit(x) |
|
| 50 |
ux <- unique(x) |
|
| 51 |
ux[which.max(tabulate(match(x, ux)))] |
|
| 52 |
} |
|
| 53 |
## aggregate to 1km resolution |
|
| 54 |
lulc2=aggregate(lulc,2,fun=function(x,na.rm=T) Mode(x)) |
|
| 55 |
## convert to factor table |
|
| 56 |
lulcf=lulc2 |
|
| 57 |
lulcf=ratify(lulcf) |
|
| 58 |
levels(lulcf) |
|
| 59 |
table(as.matrix(lulcf)) |
|
| 60 |
data(worldgrids_pal) #load palette |
|
| 61 |
IGBP=data.frame(ID=0:16,col=worldgrids_pal$IGBP[-c(18,19)], |
|
| 62 |
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)
|
|
| 63 |
IGBP$class=rownames(IGBP);rownames(IGBP)=1:nrow(IGBP) |
|
| 64 |
levels(lulcf)=list(IGBP) |
|
| 65 |
|
|
| 66 |
|
|
| 67 |
### load WORLDCLIM elevation |
|
| 68 |
dem=raster("../../tiles/h11v08/dem_h11v08.tif",format="GTiff")
|
|
| 69 |
projection(dem)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs" |
|
| 70 |
|
|
| 71 |
dif=v6-v5 |
|
| 72 |
names(dif)=month.name |
|
| 73 |
|
|
| 74 |
difm=v6m-v5m |
|
| 75 |
v5v6compare=stack(v5m,v6m,difm) |
|
| 76 |
names(v5v6compare)=c("Collection 5","Collection 6","Difference (C6-C5)")
|
|
| 77 |
|
|
| 78 |
tile=extent(v6) |
|
| 79 |
|
|
| 80 |
### compare differences between v5 and v6 by landcover type |
|
| 81 |
lulcm=as.matrix(lulc) |
|
| 82 |
forest=lulcm>=1&lulcm<=5 |
|
| 83 |
|
|
| 84 |
|
|
| 85 |
##################################### |
|
| 86 |
### compare MOD43 and MOD17 products |
|
| 87 |
|
|
| 88 |
## MOD17 |
|
| 89 |
mod17=raster("../MOD17/Npp_1km_C5.1_mean_00_to_06.tif",format="GTiff")
|
|
| 90 |
NAvalue(mod17)=32767 |
|
| 91 |
mod17=crop(projectRaster(mod17,v6,method="bilinear"),v6) |
|
| 92 |
|
|
| 93 |
mod17qc=raster("../MOD17/Npp_QC_1km_C5.1_mean_00_to_06.tif",format="GTiff")
|
|
| 94 |
mod17qc=crop(projectRaster(mod17qc,v6,method="bilinear"),v6) |
|
| 95 |
mod17qc[mod17qc<0|mod17qc>100]=NA |
|
| 96 |
|
|
| 97 |
## MOD43 via earth engine |
|
| 98 |
mod43=raster("../mod43/3b21aa90cc657523ff31e9559f36fb12.EVI_MEAN.tif",format="GTiff")
|
|
| 99 |
mod43=crop(projectRaster(mod43,v6,method="bilinear"),v6) |
|
| 100 |
|
|
| 101 |
mod43qc=raster("../mod43/3b21aa90cc657523ff31e9559f36fb12.Percent_Cloudy.tif",format="GTiff")
|
|
| 102 |
mod43qc=crop(projectRaster(mod43qc,v6,method="bilinear"),v6) |
|
| 103 |
mod43qc[mod43qc<0|mod43qc>100]=NA |
|
| 104 |
|
|
| 105 |
## Summary plot of mod17 and mod43 |
|
| 106 |
modprod=stack(mod17/cellStats(mod17,max)*100,mod17qc,mod43,mod43qc) |
|
| 107 |
names(modprod)=c("MOD17","MOD17qc","MOD43","MOD43qc")
|
|
| 108 |
|
|
| 109 |
|
|
| 110 |
### |
|
| 111 |
|
|
| 112 |
n=100 |
|
| 113 |
at=seq(0,100,len=n) |
|
| 114 |
cols=grey(seq(0,1,len=n)) |
|
| 115 |
cols=rainbow(n) |
|
| 116 |
bgyr=colorRampPalette(c("blue","green","yellow","red"))
|
|
| 117 |
cols=bgyr(n) |
|
| 118 |
|
|
| 119 |
#levelplot(lulcf,margin=F,layers="LULC") |
|
| 120 |
|
|
| 121 |
m=3 |
|
| 122 |
mcompare=stack(subset(v5,m),subset(v6,m)) |
|
| 123 |
|
|
| 124 |
mdiff=subset(v5,m)-subset(v6,m) |
|
| 125 |
names(mcompare)=c("Collection_5","Collection_6")
|
|
| 126 |
names(mdiff)=c("Collection_5-Collection_6")
|
|
| 127 |
|
|
| 128 |
|
|
| 129 |
CairoPDF("output/mod35compare.pdf",width=11,height=8.5)
|
|
| 130 |
#CairoPNG("output/mod35compare_%d.png",units="in", width=11,height=8.5,pointsize=4000,dpi=1200,antialias="subpixel")
|
|
| 131 |
|
|
| 132 |
### LANDCOVER |
|
| 133 |
levelplot(lulcf,col.regions=levels(lulcf)[[1]]$col,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) |
|
| 134 |
|
|
| 135 |
|
|
| 136 |
levelplot(mcompare,col.regions=cols,at=at,margin=F,sub="Frequency of MOD35 Clouds in March") |
|
| 137 |
#levelplot(dif,col.regions=bgyr(20),margin=F) |
|
| 138 |
levelplot(mdiff,col.regions=bgyr(100),at=seq(mdiff@data@min,mdiff@data@max,len=100),margin=F) |
|
| 139 |
|
|
| 140 |
|
|
| 141 |
boxplot(as.matrix(subset(dif,subset=1))~forest,varwidth=T,notch=T);abline(h=0) |
|
| 142 |
|
|
| 143 |
|
|
| 144 |
levelplot(modprod,main="Missing Data (%) in MOD17 (NPP) and MOD43 (BRDF Reflectance)", |
|
| 145 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at) |
|
| 146 |
|
|
| 147 |
|
|
| 148 |
|
|
| 149 |
|
|
| 150 |
levelplot(modprod,main="Missing Data (%) in MOD17 (NPP) and MOD43 (BRDF Reflectance)", |
|
| 151 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at, |
|
| 152 |
xlim=c(-7300000,-6670000),ylim=c(0,600000)) |
|
| 153 |
|
|
| 154 |
levelplot(v5m,main="Missing Data (%) in MOD17 (NPP) and MOD43 (BRDF Reflectance)", |
|
| 155 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at, |
|
| 156 |
xlim=c(-7200000,-6670000),ylim=c(0,400000),margin=F) |
|
| 157 |
|
|
| 158 |
|
|
| 159 |
levelplot(subset(v5v6compare,1:2),main="Proportion Cloudy Days (%) in Collection 5 and 6 MOD35", |
|
| 160 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at, |
|
| 161 |
margin=F) |
|
| 162 |
|
|
| 163 |
levelplot(subset(v5v6compare,1:2),main="Proportion Cloudy Days (%) in Collection 5 and 6 MOD35", |
|
| 164 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at, |
|
| 165 |
xlim=c(-7200000,-6670000),ylim=c(0,400000),margin=F) |
|
| 166 |
|
|
| 167 |
levelplot(subset(v5v6compare,1:2),main="Proportion Cloudy Days (%) in Collection 5 and 6 MOD35", |
|
| 168 |
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at, |
|
| 169 |
xlim=c(-7500000,-7200000),ylim=c(700000,1000000),margin=F) |
|
| 170 |
|
|
| 171 |
|
|
| 172 |
dev.off() |
|
| 173 |
|
|
| 174 |
### smoothing plots |
|
| 175 |
## explore smoothed version |
|
| 176 |
td=subset(v6,m) |
|
| 177 |
## build weight matrix |
|
| 178 |
s=3 |
|
| 179 |
w=matrix(1/(s*s),nrow=s,ncol=s) |
|
| 180 |
#w[s-1,s-1]=4/12; w |
|
| 181 |
td2=focal(td,w=w) |
|
| 182 |
td3=stack(td,td2) |
|
| 183 |
|
|
| 184 |
levelplot(td3,col.regions=cols,at=at,margin=F) |
|
| 185 |
|
|
| 186 |
dev.off() |
|
| 187 |
plot(stack(difm,lulc)) |
|
| 188 |
|
|
| 189 |
### ROI |
|
| 190 |
tile_ll=projectExtent(v6, "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") |
|
| 191 |
|
|
| 192 |
62,59 |
|
| 193 |
0,3 |
|
| 194 |
|
|
| 195 |
|
|
| 196 |
|
|
| 197 |
#### export KML timeseries |
|
| 198 |
library(plotKML) |
|
| 199 |
tile="h11v08" |
|
| 200 |
file=paste("summary/MOD35_",tile,".nc",sep="")
|
|
| 201 |
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=""))
|
|
| 202 |
|
|
| 203 |
v6sp=brick(paste("MOD35_",tile,".tif",sep=""))
|
|
| 204 |
v6sp=readAll(v6sp) |
|
| 205 |
|
|
| 206 |
## wasn't working with line below, perhaps Z should just be text? not date? |
|
| 207 |
v6sp=setZ(v6sp,as.Date(paste("2011-",1:12,"-15",sep="")))
|
|
| 208 |
names(v6sp)=month.name |
|
| 209 |
|
|
| 210 |
kml_open("output/mod35.kml")
|
|
| 211 |
|
|
| 212 |
|
|
| 213 |
kml_layer.RasterBrick(v6sp, |
|
| 214 |
plot.legend = TRUE, dtime = "", tz = "GMT", |
|
| 215 |
z.lim = c(0,100),colour_scale = get("colour_scale_numeric", envir = plotKML.opts))
|
|
| 216 |
# home_url = get("home_url", envir = plotKML.opts),
|
|
| 217 |
# metadata = NULL, html.table = NULL, |
|
| 218 |
# altitudeMode = "clampToGround", balloon = FALSE, |
|
| 219 |
) |
|
| 220 |
|
|
| 221 |
logo = "http://static.tumblr.com/t0afs9f/KWTm94tpm/yale_logo.png" |
|
| 222 |
kml_screen(image.file = logo, position = "UL", sname = "YALE logo",size=c(.1,.1)) |
|
| 223 |
kml_close("mod35.kml")
|
|
| 224 |
kml_compress("mod35.kml",files=c(paste(month.name,".png",sep=""),"obj_legend.png"),zip="/usr/bin/zip")
|
|
| climate/procedures/MOD35_ExtractProcessPath.r | ||
|---|---|---|
| 5 | 5 |
library(multicore) |
| 6 | 6 |
library(raster) |
| 7 | 7 |
library(spgrass6) |
| 8 |
|
|
| 8 |
library(rgeos) |
|
| 9 | 9 |
##download 3 days of modis swath data: |
| 10 | 10 |
|
| 11 | 11 |
url="ftp://ladsweb.nascom.nasa.gov/allData/51/MOD35_L2/2012/" |
| ... | ... | |
| 18 | 18 |
## get MODLAND tile to serve as base |
| 19 | 19 |
#system("wget http://e4ftl01.cr.usgs.gov/MOLT/MOD13A3.005/2000.02.01/MOD13A3.A2000032.h00v08.005.2006271174446.hdf")
|
| 20 | 20 |
#t=raster(paste("HDF4_EOS:EOS_GRID:\"",getwd(),"/MOD13A3.A2000032.h00v08.005.2006271174446.hdf\":MOD_Grid_monthly_1km_VI:1 km monthly NDVI",sep=""))
|
| 21 |
t=raster(paste("../MOD17/Npp_1km_C5.1_mean_00_to_06.tif",sep=""))
|
|
| 21 |
t=raster(paste("../MOD17/MOD17A3_Science_NPP_mean_00_12.tif",sep=""))
|
|
| 22 |
projection(t) |
|
| 22 | 23 |
|
| 23 | 24 |
## make global extent |
| 24 | 25 |
pmodis="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs" |
| 25 | 26 |
|
| 26 | 27 |
glb=t |
| 27 |
values(glb)=NA |
|
| 28 |
#values(glb)=NA
|
|
| 28 | 29 |
glb=extend(glb,extent(-180,180,-90,90)) |
| 29 | 30 |
|
| 30 | 31 |
#glb=raster(glb,crs="+proj=longlat +datum=WGS84",nrows=42500,ncols=85000) |
| ... | ... | |
| 38 | 39 |
|
| 39 | 40 |
#### Grid and mosaic the swath data |
| 40 | 41 |
|
| 41 |
stitch="sudo MRTDATADIR=\"/usr/local/heg/data\" PGSHOME=/usr/local/heg/TOOLKIT_MTD PWD=/home/adamw /usr/local/heg/bin/swtif"
|
|
| 42 |
stitch="sudo MRTDATADIR=\"/usr/local/heg/data\" PGSHOME=/usr/local/heg/TOOLKIT_MTD PWD=/home/adamw /usr/local/heg/bin/swtif" |
|
| 42 | 43 |
|
| 43 | 44 |
files=paste(getwd(),"/",list.files("swath",pattern="hdf$",full=T),sep="")
|
| 44 | 45 |
|
| 45 | 46 |
## vars to process |
| 46 | 47 |
vars=as.data.frame(matrix(c( |
| 47 | 48 |
"Cloud_Mask", "CM", |
| 48 |
# "Quality_Assurance", "QA"), |
|
| 49 |
byrow=T,ncol=2,dimnames=list(1:2,c("variable","varid"))),stringsAsFactors=F)
|
|
| 50 |
|
|
| 51 |
## establish sudo priveleges to run swtif |
|
| 52 |
system("sudo ls")
|
|
| 53 |
|
|
| 54 |
mclapply(files,function(file){
|
|
| 55 |
|
|
| 56 |
tempfile=paste(tempdir(),"/",basename(file),sep="") |
|
| 57 |
# tempfile2=paste("gridded/",sub("hdf","tif",basename(tempfile)),sep="")
|
|
| 58 |
tempfile2=paste(tempdir(),"/",sub("hdf","tif",basename(tempfile)),sep="")
|
|
| 59 |
|
|
| 60 |
outfile=paste("gridded2/",sub("hdf","tif",basename(tempfile)),sep="")
|
|
| 61 |
|
|
| 62 |
if(file.exists(outfile)) return(c(file,0)) |
|
| 63 |
## First write the parameter file (careful, heg is very finicky!) |
|
| 64 |
hdr=paste("NUM_RUNS = ",length(vars$varid),"|MULTI_BAND_HDFEOS:",length(vars$varid),sep="")
|
|
| 65 |
grp=paste("
|
|
| 49 |
"Sensor_Azimuth", "ZA", |
|
| 50 |
"Sensor_Zenith", "SZ"), |
|
| 51 |
byrow=T,ncol=2,dimnames=list(1:3,c("variable","varid"))),stringsAsFactors=F)
|
|
| 52 |
|
|
| 53 |
## global bounding box |
|
| 54 |
gbb=cbind(c(-180,-180,180,180,-180),c(-85,85,85,-85,-85)) |
|
| 55 |
gpp = SpatialPolygons(list(Polygons(list(Polygon(gbb)),1))) |
|
| 56 |
proj4string(gpp)=projection(glb) |
|
| 57 |
|
|
| 58 |
|
|
| 59 |
getpath<- function(file){
|
|
| 60 |
setwd(tempdir()) |
|
| 61 |
bfile=sub(".hdf","",basename(file))
|
|
| 62 |
tempfile_path=paste(tempdir(),"/path_",basename(file),sep="") #gridded path |
|
| 63 |
tempfile_sz=paste(tempdir(),"/sz_",basename(file),sep="") # gridded sensor zenith |
|
| 64 |
tempfile2_path=paste(tempdir(),"/",bfile,".tif",sep="") #gridded/masked/processed path |
|
| 65 |
outfile=paste("~/acrobates/adamw/projects/interp/data/modis/mod35/gridded/",bfile,".tif",sep="") #final file
|
|
| 66 |
if(file.exists(outfile)) return(c(file,0)) |
|
| 67 |
## get bounding coordinates |
|
| 68 |
glat=as.numeric(do.call(c,strsplit(sub("GRINGPOINTLATITUDE=","",system(paste("gdalinfo ",file," | grep GRINGPOINTLATITUDE"),intern=T)),split=",")))
|
|
| 69 |
glon=as.numeric(do.call(c,strsplit(sub("GRINGPOINTLONGITUDE=","",system(paste("gdalinfo ",file," | grep GRINGPOINTLONGITUDE"),intern=T)),split=",")))
|
|
| 70 |
bb=cbind(c(glon,glon[1]),c(glat,glat[1])) |
|
| 71 |
pp = SpatialPolygons(list(Polygons(list(Polygon(bb)),1))) |
|
| 72 |
proj4string(pp)=projection(glb) |
|
| 73 |
ppc=gIntersection(pp,gpp) |
|
| 74 |
ppc=gBuffer(ppc,width=0.3) #buffer a little to remove gaps between images |
|
| 75 |
## First write the parameter file (careful, heg is very finicky!) |
|
| 76 |
hdr=paste("NUM_RUNS = ",length(vars$varid),"|MULTI_BAND_HDFEOS:",length(vars$varid),sep="")
|
|
| 77 |
grp=paste("
|
|
| 66 | 78 |
BEGIN |
| 67 | 79 |
INPUT_FILENAME=",file," |
| 68 | 80 |
OBJECT_NAME=mod35 |
| ... | ... | |
| 70 | 82 |
BAND_NUMBER = ",1:length(vars$varid)," |
| 71 | 83 |
OUTPUT_PIXEL_SIZE_X=0.008333333 |
| 72 | 84 |
OUTPUT_PIXEL_SIZE_Y=0.008333333 |
| 73 |
SPATIAL_SUBSET_UL_CORNER = ( 90 -180 )
|
|
| 74 |
SPATIAL_SUBSET_LR_CORNER = ( -90 180 )
|
|
| 85 |
SPATIAL_SUBSET_UL_CORNER = ( ",bbox(ppc)[2,2]," ",bbox(ppc)[1,1]," )
|
|
| 86 |
SPATIAL_SUBSET_LR_CORNER = ( ",bbox(ppc)[2,1]," ",bbox(ppc)[1,2]," )
|
|
| 75 | 87 |
OUTPUT_OBJECT_NAME = mod35| |
| 76 | 88 |
RESAMPLING_TYPE =NN |
| 77 | 89 |
OUTPUT_PROJECTION_TYPE = GEO |
| ... | ... | |
| 80 | 92 |
# projection parameters from http://landweb.nascom.nasa.gov/cgi-bin/QA_WWW/newPage.cgi?fileName=sn_gctp |
| 81 | 93 |
ELLIPSOID_CODE = WGS84 |
| 82 | 94 |
OUTPUT_TYPE = HDFEOS |
| 83 |
OUTPUT_FILENAME= ",tempfile," |
|
| 95 |
OUTPUT_FILENAME= ",tempfile_path,"
|
|
| 84 | 96 |
END |
| 85 | 97 |
|
| 86 | 98 |
",sep="") |
| 87 | 99 |
|
| 88 |
## if any remnants from previous runs remain, delete them |
|
| 89 |
if(length(list.files(tempdir(),pattern=basename(file)))>0) |
|
| 90 |
file.remove(list.files(tempdir(),pattern=basename(file),full=T)) |
|
| 91 |
## write it to a file |
|
| 92 |
cat( c(hdr,grp) , file=paste(tempdir(),"/",basename(file),"_MODparms.txt",sep="")) |
|
| 93 |
|
|
| 94 |
## now run the swath2grid tool |
|
| 95 |
## write the gridded file |
|
| 96 |
print(paste("Starting",file))
|
|
| 97 |
system(paste("(",stitch," -p ",tempdir(),"/",basename(file),"_MODparms.txt -d ; echo $$)",sep=""),intern=T,ignore.stderr=F)
|
|
| 98 |
## convert to land path |
|
| 99 |
d=raster(paste("HDF4_EOS:EOS_GRID:\"",tempfile,"\":mod35:Cloud_Mask_0",sep=""))
|
|
| 100 |
# za=raster(paste("HDF4_EOS:EOS_GRID:\"",tempfile,"\":mod35:Quality_Assurance_1",sep=""))
|
|
| 101 |
## add projection information - ellipsoid should be wgs84 |
|
| 102 |
projection(d)=projection(glb) |
|
| 103 |
extent(d)=alignExtent(d,extent(glb)) |
|
| 104 |
# d=readAll(d) |
|
| 105 |
getlc=function(x) {(x%/%2^6) %% 2^2}
|
|
| 106 |
calc(d,fun=getlc,filename=outfile,options=c("COMPRESS=LZW", "LEVEL=9","PREDICTOR=2"),datatype="INT1U",overwrite=T)
|
|
| 107 |
|
|
| 108 |
### warp it to align all pixels |
|
| 109 |
system(paste("gdalwarp -overwrite -tap -tr 0.008333333 0.008333333 -co COMPRESS=LZW -co LEVEL=9 -co PREDICTOR=2 -s_srs \"",projection(glb),"\" ",tempfile2," ",outfile,sep=""))
|
|
| 110 |
|
|
| 111 |
# getqa=function(x) {(x%/%2^1) %% 2^3}
|
|
| 112 |
# za=calc(za,fun=getqa)#,filename=outfile,options=c("COMPRESS=LZW", "LEVEL=9","PREDICTOR=2"),datatype="INT1U",overwrite=T)
|
|
| 113 |
## resample to line up with glb so we can use mosaic later |
|
| 114 |
## delete temporary files |
|
| 115 |
file.remove(tempfile,tempfile2) |
|
| 116 |
return(c(file,1)) |
|
| 117 |
}) |
|
| 100 |
## if any remnants from previous runs remain, delete them |
|
| 101 |
if(length(list.files(tempdir(),pattern=bfile)>0)) |
|
| 102 |
file.remove(list.files(tempdir(),pattern=bfile,full=T)) |
|
| 103 |
## write it to a file |
|
| 104 |
cat( c(hdr,grp) , file=paste(tempdir(),"/",basename(file),"_MODparms.txt",sep="")) |
|
| 105 |
## now run the swath2grid tool |
|
| 106 |
## write the gridded file |
|
| 107 |
print(paste("Starting",file))
|
|
| 108 |
system(paste("(",stitch," -p ",tempdir(),"/",basename(file),"_MODparms.txt -d)",sep=""),intern=F)#,ignore.stderr=F)
|
|
| 109 |
############## Now run the 5km summary |
|
| 110 |
## First write the parameter file (careful, heg is very finicky!) |
|
| 111 |
hdr=paste("NUM_RUNS = ",nrow(vars[-1,]),"|MULTI_BAND_HDFEOS:",nrow(vars[-1,]),sep="")
|
|
| 112 |
grp=paste("
|
|
| 113 |
BEGIN |
|
| 114 |
INPUT_FILENAME=",file," |
|
| 115 |
OBJECT_NAME=mod35 |
|
| 116 |
FIELD_NAME=",vars$variable[-1],"| |
|
| 117 |
BAND_NUMBER = ",1," |
|
| 118 |
OUTPUT_PIXEL_SIZE_X=0.008333333 |
|
| 119 |
#0.0416666 |
|
| 120 |
OUTPUT_PIXEL_SIZE_Y=0.008333333 |
|
| 121 |
#0.0416666 |
|
| 122 |
SPATIAL_SUBSET_UL_CORNER = ( ",bbox(ppc)[2,2]," ",bbox(ppc)[1,1]," ) |
|
| 123 |
SPATIAL_SUBSET_LR_CORNER = ( ",bbox(ppc)[2,1]," ",bbox(ppc)[1,2]," ) |
|
| 124 |
#OUTPUT_OBJECT_NAME = mod35| |
|
| 125 |
RESAMPLING_TYPE =NN |
|
| 126 |
OUTPUT_PROJECTION_TYPE = GEO |
|
| 127 |
OUTPUT_PROJECTION_PARAMETERS = ( 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ) |
|
| 128 |
#OUTPUT_PROJECTION_PARAMETERS = ( 6371007.181 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ) |
|
| 129 |
# projection parameters from http://landweb.nascom.nasa.gov/cgi-bin/QA_WWW/newPage.cgi?fileName=sn_gctp |
|
| 130 |
ELLIPSOID_CODE = WGS84 |
|
| 131 |
OUTPUT_TYPE = HDFEOS |
|
| 132 |
OUTPUT_FILENAME= ",tempfile_sz," |
|
| 133 |
END |
|
| 134 |
|
|
| 135 |
",sep="") |
|
| 136 |
|
|
| 137 |
## write it to a file |
|
| 138 |
cat( c(hdr,grp) , file=paste(tempdir(),"/",basename(file),"_MODparms_angle.txt",sep="")) |
|
| 139 |
|
|
| 140 |
## now run the swath2grid tool |
|
| 141 |
## write the gridded file |
|
| 142 |
print(paste("Starting",file))
|
|
| 143 |
system(paste("(",stitch," -p ",tempdir(),"/",basename(file),"_MODparms_angle.txt -d)",sep=""),intern=F,ignore.stderr=F)
|
|
| 144 |
####### import to R for processing |
|
| 145 |
if(!file.exists(tempfile_path)) {
|
|
| 146 |
file.remove(list.files(tempdir(),pattern=bfile,full=T)) |
|
| 147 |
return(c(file,0)) |
|
| 148 |
} |
|
| 149 |
## convert to land path |
|
| 150 |
d=raster(paste("HDF4_EOS:EOS_GRID:\"",tempfile_path,"\":mod35:Cloud_Mask_0",sep=""))
|
|
| 151 |
sz=raster(paste("HDF4_EOS:EOS_GRID:\"",tempfile_sz,"\":mod35:Sensor_Zenith",sep=""))
|
|
| 152 |
NAvalue(sz)=0 |
|
| 153 |
## resample sensor angles to 1km grid and mask paths with angles >=30 |
|
| 154 |
# sz2=resample(sz,d,method="ngb",file=sub("hdf","tif",tempfile_sz),overwrite=T)
|
|
| 155 |
getlc=function(x,y) {ifelse(y==0|y>40,NA,((x%/%2^6) %% 2^2))}
|
|
| 156 |
path= overlay(d,sz,fun=getlc,filename=tempfile2_path,options=c("COMPRESS=LZW", "LEVEL=9","PREDICTOR=2"),datatype="INT1U",overwrite=T)
|
|
| 157 |
### warp them to align all pixels |
|
| 158 |
system(paste("gdalwarp -overwrite -srcnodata 255 -dstnodata 255 -tap -tr 0.008333333 0.008333333 -co COMPRESS=LZW -co ZLEVEL=9 -co PREDICTOR=2 -s_srs \"",projection(t),"\" ",tempfile2_path," ",outfile,sep=""))
|
|
| 159 |
## delete temporary files |
|
| 160 |
file.remove(list.files(tempdir(),pattern=bfile,full=T)) |
|
| 161 |
return(c(file,1)) |
|
| 162 |
} |
|
| 163 |
|
|
| 164 |
|
|
| 165 |
## establish sudo priveleges to run swtif |
|
| 166 |
system("sudo ls"); mclapply(files,getpath,mc.cores=10)
|
|
| 118 | 167 |
|
| 119 | 168 |
## check gdal can read all of them |
| 120 | 169 |
gfiles=list.files("gridded",pattern="tif$",full=T)
|
| ... | ... | |
| 130 | 179 |
|
| 131 | 180 |
file.remove(gfiles[check==0]) |
| 132 | 181 |
|
| 133 |
# origin(raster(gfiles[5])) |
|
| 134 |
|
|
| 135 |
|
|
| 136 |
system(paste("gdalinfo ",gfiles[1]," | grep Origin"))
|
|
| 137 |
system(paste("gdalinfo ",gfiles[2]," | grep Origin"))
|
|
| 138 |
system(paste("gdalinfo ",gfiles[3]," | grep Origin"))
|
|
| 139 |
|
|
| 140 |
system(paste("gdalwarp -tap -tr 0.008333333 0.008333333 -s_srs \"",projection(glb),"\" ",gfiles[1]," test1.tif"))
|
|
| 141 |
system(paste("gdalwarp -tap -tr 0.008333333 0.008333333 -s_srs \"",projection(glb),"\" ",gfiles[2]," test2.tif"))
|
|
| 142 |
system(paste("gdalinfo test1.tif | grep Origin"))
|
|
| 143 |
system(paste("gdalinfo test2.tif | grep Origin"))
|
|
| 144 |
system(paste("pkmosaic ",paste("-i",gfiles[1:2],collapse=" ")," -o MOD35_path2.tif -m 6 -v -t 255"))
|
|
| 145 |
system(paste("pkmosaic ",paste("-i",c("test1.tif","test2.tif"),collapse=" ")," -o MOD35_path2.tif -m 6 -v -max 10 -ot Byte"))
|
|
| 182 |
## use new gdal |
|
| 183 |
system(paste("/usr/local/gdal-1.10.0/bin/gdalwarp -wm 900 -overwrite -co COMPRESS=LZW -co PREDICTOR=2 -multi -r mode gridded/*.tif MOD35_path_gdalwarp.tif"))
|
|
| 146 | 184 |
|
| 147 | 185 |
|
| 186 |
# origin(raster(gfiles[5])) |
|
| 187 |
|
|
| 148 | 188 |
## try with pktools |
| 149 | 189 |
## global |
| 150 |
system(paste("pkmosaic -co COMPRESS=LZW -co PREDICTOR=2 ",paste("-i",list.files("gridded",full=T,pattern="tif$"),collapse=" ")," -o MOD35_path_pkmosaic.tif -m 6 -v -t 255 -t 0 &"))
|
|
| 190 |
system(paste("pkmosaic -co COMPRESS=LZW -co PREDICTOR=2 ",paste("-i",list.files("gridded",full=T,pattern="tif$"),collapse=" ")," -o MOD35_path_pkmosaic_max.tif -m 2 -v -t 255 -t 0 &"))
|
|
| 151 | 191 |
#bb="-ulx -180 -uly 90 -lrx 180 -lry -90" |
| 152 | 192 |
#bb="-ulx -180 -uly 90 -lrx 170 -lry 80" |
| 153 | 193 |
bb="-ulx -72 -uly 11 -lrx -59 -lry -1" |
| 154 | 194 |
|
| 155 | 195 |
|
| 156 | 196 |
#expand.grid(x=seq(-180,170,by=10),y=seq(-90,80)) |
| 157 |
|
|
| 158 |
system(paste("pkmosaic ",bb," -co COMPRESS=LZW -co PREDICTOR=2 ",paste("-i",list.files("gridded",full=T,pattern="tif$"),collapse=" ")," -o h11v08_path_pkmosaic.tif -m 6 -v -t 255 -t 0 &"))
|
|
| 197 |
gf2= grep("2012009[.]03",gfiles,value=T)
|
|
| 198 |
system(paste("pkmosaic ",bb," -co COMPRESS=LZW -co PREDICTOR=2 ",paste("-i",gf2,collapse=" ")," -o h11v08_path_pkmosaic.tif -ot Byte -m 7 -v -t 255"))
|
|
| 159 | 199 |
|
| 160 | 200 |
# bounding box? |
| 161 | 201 |
|
| ... | ... | |
| 193 | 233 |
unlink_.gislock() |
| 194 | 234 |
system(paste("rm -frR ",tf,sep=""))
|
| 195 | 235 |
|
| 196 |
|
|
| 197 | 236 |
######################### |
| 198 | 237 |
|
| 199 | 238 |
|
| climate/procedures/MODIS_Cloud.py | ||
|---|---|---|
| 1 |
import ee |
|
| 2 |
from ee import mapclient |
|
| 3 |
|
|
| 4 |
import logging |
|
| 5 |
logging.basicConfig() |
|
| 6 |
|
|
| 7 |
MY_SERVICE_ACCOUNT = '511722844190@developer.gserviceaccount.com' # replace with your service account |
|
| 8 |
MY_PRIVATE_KEY_FILE = '/Users/adamw/EarthEngine-privatekey.p12' # replace with you private key file path |
|
| 9 |
|
|
| 10 |
ee.Initialize(ee.ServiceAccountCredentials(MY_SERVICE_ACCOUNT, MY_PRIVATE_KEY_FILE)) |
|
| 11 |
|
|
| 12 |
#// EVI_Cloud_month |
|
| 13 |
|
|
| 14 |
#// Make land mask |
|
| 15 |
#// Select the forest classes from the MODIS land cover image and intersect them |
|
| 16 |
#// with elevations above 1000m. |
|
| 17 |
elev = ee.Image('srtm90_v4');
|
|
| 18 |
cover = ee.Image('MCD12Q1/MCD12Q1_005_2001_01_01').select('Land_Cover_Type_1');
|
|
| 19 |
blank = ee.Image(0); |
|
| 20 |
#// Where (cover == 0) and (elev > 0), set the output to 1. |
|
| 21 |
land = blank.where( |
|
| 22 |
cover.neq(0).and(cover.neq(15)),//.and(elev.gt(0)), |
|
| 23 |
1); |
|
| 24 |
|
|
| 25 |
palette = ["aec3d4", // water |
|
| 26 |
"152106", "225129", "369b47", "30eb5b", "387242", // forest |
|
| 27 |
"6a2325", "c3aa69", "b76031", "d9903d", "91af40", // shrub, grass, savanah |
|
| 28 |
"111149", // wetlands |
|
| 29 |
"cdb33b", // croplands |
|
| 30 |
"cc0013", // urban |
|
| 31 |
"33280d", // crop mosaic |
|
| 32 |
"d7cdcc", // snow and ice |
|
| 33 |
"f7e084", // barren |
|
| 34 |
"6f6f6f"].join(',');// tundra
|
|
| 35 |
|
|
| 36 |
#// make binary map for forest/nonforest |
|
| 37 |
forest = blank.where(cover.gte(1).and(cover.lte(5)),1); |
|
| 38 |
|
|
| 39 |
#//addToMap(forest, {min: 0, max: 1});
|
|
| 40 |
#//addToMap(cover, {min: 0, max: 17, palette:palette});
|
|
| 41 |
|
|
| 42 |
#// MODIS EVI Collection |
|
| 43 |
collection = ee.ImageCollection("MCD43A4_EVI");
|
|
| 44 |
|
|
| 45 |
#//define reducers and filters |
|
| 46 |
COUNT = ee.call("Reducer.count");
|
|
| 47 |
#// Loop through months and get monthly % cloudiness |
|
| 48 |
#//for (var month = 1; month < 2; month += 1) {
|
|
| 49 |
#//month=1; |
|
| 50 |
FilterMonth = ee.Filter(ee.call("Filter.calendarRange",
|
|
| 51 |
start=month,end=month,field="month")); |
|
| 52 |
tmonth=collection.filter(FilterMonth) |
|
| 53 |
|
|
| 54 |
n=tmonth.getInfo().features.length; #// Get total number of images |
|
| 55 |
print(n+" Layers in the collection for month "+month) |
|
| 56 |
tcloud=tmonth.reduce(COUNT).float() |
|
| 57 |
c=ee.Image(n); #// make raster with constant value of n |
|
| 58 |
c1=ee.Image(-1); #// make raster with constant value of -1 to convert to % cloudy |
|
| 59 |
|
|
| 60 |
#///////////////////////////////////////////////// |
|
| 61 |
#// Generate the cloud frequency surface: |
|
| 62 |
#// 1 Calculate the number of days with measurements |
|
| 63 |
#// 2 divide by the total number of layers |
|
| 64 |
#// 3 Add -1 and multiply by -1 to invert to % cloudy |
|
| 65 |
#// 4 Rename to "PCloudy_month" |
|
| 66 |
tcloud = tcloud.divide(c).add(c1).multiply(c1).expression("b()*100").int8().select_([0],["PCloudy_"+month]);
|
|
| 67 |
#//if(month==1) {var cloud=tcloud} // if first year, make new object
|
|
| 68 |
#//if(month>1) {var cloud=cloud.addBands(tcloud)} // if not first year, append to first year
|
|
| 69 |
#//} //end loop over months |
|
| 70 |
# // end loop over years |
|
| 71 |
|
|
| 72 |
|
|
| 73 |
#//print(evi_miss.stats()) |
|
| 74 |
#//addToMap(tcloud,{min:0,max:100,palette:"000000,00FF00,FF0000"});
|
|
| 75 |
#//addToMap(elev,{min:0,max:2500,palette:"000000,00FF00,FF0000"});
|
|
| 76 |
|
|
| 77 |
#//centerMap(-122.418, 37.72, 10); |
|
| 78 |
path = tcloud.getDownloadURL({
|
|
| 79 |
'scale': 1000, |
|
| 80 |
'crs': 'EPSG:4326', |
|
| 81 |
'region': '[[-90, 0], [-90, 20], [-50, 0], [-50, 20]]' //h11v08 |
|
| 82 |
#// 'region': '[[-180, -90], [-180, 90], [180, -90], [180, 90]]' |
|
| 83 |
}); |
|
| 84 |
print('https://earthengine.googleapis.com' + path);
|
|
| 85 |
|
|
| climate/procedures/NDP-026D.R | ||
|---|---|---|
| 5 | 5 |
library(multicore) |
| 6 | 6 |
library(latticeExtra) |
| 7 | 7 |
library(doMC) |
| 8 |
library(raster) |
|
| 8 |
library(rasterVis)
|
|
| 9 | 9 |
library(rgdal) |
| 10 | 10 |
## register parallel processing |
| 11 | 11 |
registerDoMC(20) |
| ... | ... | |
| 105 | 105 |
coordinates(cldms)=c("lon","lat")
|
| 106 | 106 |
projection(cldms)=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
|
| 107 | 107 |
|
| 108 |
##make spatial object |
|
| 109 |
cldys=cldy |
|
| 110 |
coordinates(cldys)=c("lon","lat")
|
|
| 111 |
projection(cldys)=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
|
|
| 112 |
|
|
| 108 | 113 |
#### Evaluate MOD35 Cloud data |
| 109 | 114 |
mod35=brick("../modis/mod35/MOD35_h11v08.nc",varname="CLD01")
|
| 110 | 115 |
mod35sd=brick("../modis/mod35/MOD35_h11v08.nc",varname="CLD_sd")
|
| ... | ... | |
| 112 | 117 |
|
| 113 | 118 |
|
| 114 | 119 |
### use data from google earth engine |
| 115 |
mod35=brick("../modis/mod09/global_2009/")
|
|
| 120 |
mod35=raster("../modis/mod09/global_2009/MOD35_2009.tif")
|
|
| 116 | 121 |
mod09=raster("../modis/mod09/global_2009/MOD09_2009.tif")
|
| 117 | 122 |
|
| 123 |
## LULC |
|
| 124 |
system(paste("gdalwarp -r near -co \"COMPRESS=LZW\" -tr ",paste(res(mod09),collapse=" ",sep=""),
|
|
| 125 |
"-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")) |
|
| 126 |
lulc=raster("../modis/mod12/MCD12Q1_IGBP_2005_v51.tif")
|
|
| 127 |
lulc=ratify(lulc) |
|
| 128 |
require(plotKML); data(worldgrids_pal) #load IGBP palette |
|
| 129 |
IGBP=data.frame(ID=0:16,col=worldgrids_pal$IGBP[-c(18,19)],stringsAsFactors=F) |
|
| 130 |
IGBP$class=rownames(IGBP);rownames(IGBP)=1:nrow(IGBP) |
|
| 131 |
levels(lulc)=list(IGBP) |
|
| 132 |
lulc=crop(lulc,mod09) |
|
| 133 |
|
|
| 118 | 134 |
n=100 |
| 119 | 135 |
at=seq(0,100,length=n) |
| 120 | 136 |
colr=colorRampPalette(c("black","green","red"))
|
| 121 | 137 |
cols=colr(n) |
| 122 | 138 |
|
| 123 |
levelplot(mod09,col.regions=cols,at=at) |
|
| 139 |
dif=mod35-mod09 |
|
| 140 |
bwplot(dif~as.factor(lulc)) |
|
| 124 | 141 |
|
| 142 |
levelplot(mod35,col.regions=cols,at=at,margins=F,maxpixels=1e6)#,xlim=c(-100,-50),ylim=c(0,10)) |
|
| 143 |
levelplot(lulc,att="class",col.regions=levels(lulc)[[1]]$col,margin=F,maxpixels=1e6) |
|
| 125 | 144 |
|
| 126 |
cldms=spTransform(cldms,CRS(projection(mod35)))
|
|
| 145 |
#cldys=spTransform(cldys,CRS(projection(mod35)))
|
|
| 127 | 146 |
|
| 128 | 147 |
mod35v=foreach(m=unique(cldm$month),.combine="rbind") %do% {
|
| 129 | 148 |
dr=subset(mod35,subset=m);projection(dr)=projection(mod35) |
| ... | ... | |
| 134 | 153 |
print(m) |
| 135 | 154 |
return(ds@data[!is.na(ds$mod35),])} |
| 136 | 155 |
|
| 156 |
y=2009 |
|
| 157 |
d=cldys[cldys$year==y,] |
|
| 158 |
|
|
| 159 |
d$mod35_10=unlist(extract(mod35,d,buffer=10000,fun=mean,na.rm=T)) |
|
| 160 |
d$mod09_10=unlist(extract(mod09,d,buffer=10000,fun=mean,na.rm=T)) |
|
| 161 |
d$dif=d$mod35_10-d$mod09_10 |
|
| 162 |
d$dif2=d$mod35_10-d$cld |
|
| 163 |
|
|
| 164 |
d$lulc=unlist(extract(lulc,d)) |
|
| 165 |
d$lulc_10=unlist(extract(lulc,d,buffer=10000,fun=mode,na.rm=T)) |
|
| 166 |
d$lulc=factor(d$lulc,labels=IGBP$class) |
|
| 167 |
|
|
| 168 |
save(d,file="annualsummary.Rdata") |
|
| 169 |
|
|
| 170 |
## quick model to explore fit |
|
| 171 |
plot(cld~mod35,groups=lulc,data=d) |
|
| 172 |
summary(lm(cld~mod35+as.factor(lulc),data=d)) |
|
| 173 |
summary(lm(cld~mod09_10,data=d)) |
|
| 174 |
summary(lm(cld~mod09_10+as.factor(lulc),data=d)) |
|
| 175 |
summary(lm(cld~mod09_10+as.factor(lulc),data=d)) |
|
| 176 |
|
|
| 177 |
### exploratory plots |
|
| 178 |
xyplot(cld~mod09_10,groups=lulc,data=d@data,pch=16,cex=.5)+layer(panel.abline(0,1,col="red")) |
|
| 179 |
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"))
|
|
| 180 |
xyplot(cld~mod35_10|as.factor(lulc),data=d@data,pch=16,cex=.5)+layer(panel.abline(0,1,col="red")) |
|
| 181 |
xyplot(mod35_10~mod09_10|as.factor(lulc),data=d@data,pch=16,cex=.5)+layer(panel.abline(0,1,col="red")) |
|
| 182 |
|
|
| 183 |
densityplot(stack(mod35,mod09)) |
|
| 184 |
boxplot(mod35,lulc) |
|
| 185 |
|
|
| 186 |
bwplot(mod09~mod35|cut(y,5),data=stack(mod09,mod35)) |
|
| 187 |
|
|
| 137 | 188 |
## month factors |
| 138 | 189 |
cldm$month2=factor(cldm$month,labels=month.name) |
| 139 | 190 |
## add a color key |
Also available in: Unified diff
Updated script to process MOD35 Processing Path using HEG to interpolate sensor zenith observations. Also working on figures for MOD35 C5 vs C6 comparison in MOD35C5_Evaluation script