Revision a57c4e3d
Added by Adam Wilson over 11 years ago
| climate/procedures/MOD35C5_Evaluation.r | ||
|---|---|---|
| 10 | 10 |
library(reshape) |
| 11 | 11 |
library(rgeos) |
| 12 | 12 |
|
| 13 |
## landcover |
|
| 14 |
if(!file.exists("data/MCD12Q1_IGBP_2005_v51_1km_wgs84.tif")){
|
|
| 15 |
system(paste("gdalwarp -r near -ot Byte -co \"COMPRESS=LZW\"",
|
|
| 16 |
" ~/acrobatesroot/jetzlab/Data/environ/global/landcover/MODIS/MCD12Q1_IGBP_2005_v51.tif ", |
|
| 17 |
" -t_srs \"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs\" ", |
|
| 18 |
tempdir(),"/MCD12Q1_IGBP_2005_v51_wgs84.tif -overwrite ",sep="")) |
|
| 19 |
lulc=raster(paste(tempdir(),"/MCD12Q1_IGBP_2005_v51_wgs84.tif",sep="")) |
|
| 20 |
## aggregate to 1km resolution |
|
| 21 |
lulc2=aggregate(lulc,2,fun=function(x,na.rm=T) {
|
|
| 22 |
x=na.omit(x) |
|
| 23 |
ux <- unique(x) |
|
| 24 |
ux[which.max(tabulate(match(x, ux)))] |
|
| 25 |
},file=paste(tempdir(),"/1km.tif",sep="")) |
|
| 26 |
writeRaster(lulc2,"data/MCD12Q1_IGBP_2005_v51_1km_wgs84.tif",options=c("COMPRESS=LZW","ZLEVEL=9","PREDICTOR=2"),datatype="INT1U",overwrite=T)
|
|
| 27 |
} |
|
| 28 |
## read it in |
|
| 29 |
lulc=raster("data/MCD12Q1_IGBP_2005_v51_1km_wgs84.tif")
|
|
| 30 |
# lulc=ratify(lulc) |
|
| 31 |
data(worldgrids_pal) #load palette |
|
| 32 |
IGBP=data.frame(ID=0:16,col=worldgrids_pal$IGBP[-c(18,19)], |
|
| 33 |
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)
|
|
| 34 |
IGBP$class=rownames(IGBP);rownames(IGBP)=1:nrow(IGBP) |
|
| 35 |
levels(lulc)=list(IGBP) |
|
| 36 |
extent(lulc)=alignExtent(lulc,mod09) |
|
| 37 |
names(lulc)="MCD12Q1" |
|
| 38 |
|
|
| 39 |
## make land mask |
|
| 40 |
if(!file.exists("data/land.tif"))
|
|
| 41 |
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)
|
|
| 42 |
land=raster("data/land.tif")
|
|
| 43 |
|
|
| 44 | 13 |
## get % cloudy |
| 45 | 14 |
mod09=raster("data/MOD09_2009.tif")
|
| 46 | 15 |
names(mod09)="MOD09CF" |
| 47 | 16 |
NAvalue(mod09)=0 |
| 48 | 17 |
|
| 49 | 18 |
mod35c5=raster("data/MOD35_2009.tif")
|
| 50 |
mod35c5=crop(mod35c5,mod09) |
|
| 51 | 19 |
names(mod35c5)="C5MOD35CF" |
| 52 | 20 |
NAvalue(mod35c5)=0 |
| 53 | 21 |
|
| 54 |
system("gdalinfo /home/adamw/acrobates/adamw/projects/interp/data/modis/mod35/summary/.nc")
|
|
| 55 |
|
|
| 56 | 22 |
## mod35C6 annual |
| 57 |
system("/usr/local/gdal-1.10.0/bin/gdalbuildvrt -a_srs '+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs' -sd 1 -b 1 data/MOD35C6.vrt /home/adamw/acrobates/adamw/projects/interp/data/modis/mod35/summary/*2009mean.nc ")
|
|
| 58 |
system("align.sh data/MOD35C6.vrt data/MOD35_2009.tif data/MOD35C6_2009.tif")
|
|
| 23 |
if(!file.exists("data/MOD35C6_2009.tif")){
|
|
| 24 |
system("/usr/local/gdal-1.10.0/bin/gdalbuildvrt -a_srs '+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs' -sd 1 -b 1 data/MOD35C6.vrt /home/adamw/acrobates/adamw/projects/interp/data/modis/mod35/summary/*2009mean.nc ")
|
|
| 25 |
system("align.sh data/MOD35C6.vrt data/MOD09_2009.tif data/MOD35C6_2009.tif")
|
|
| 26 |
} |
|
| 59 | 27 |
mod35c6=raster("data/MOD35C6_2009.tif")
|
| 60 | 28 |
mod35c6=crop(mod35c6,mod09) |
| 61 | 29 |
names(mod35c6)="C6MOD35CF" |
| 62 | 30 |
NAvalue(mod35c6)=255 |
| 63 | 31 |
|
| 64 |
## Monthly |
|
| 65 |
#system("find /home/adamw/acrobates/adamw/projects/interp/data/modis/mod35/summary/ -name '*[0-9].nc' > data/summaryfiles.txt")
|
|
| 66 |
#system("/usr/local/gdal-1.10.0/bin/gdalbuidvrt -allow_projection_difference -b 1 -sd 1 data/MOD35C6_monthly.vrt ")
|
|
| 67 |
##system("/usr/local/gdal-1.10.0/bin/gdalwarp --optfile data/summaryfiles.txt data/MOD35C6_monthly.tif ")#
|
|
| 68 |
#system("/usr/local/gdal-1.10.0/bin/gdal_translate -stats -co \"COMPRESS=LZW\" -of GTiff data/MOD35C6_monthly.vrt data/MOD35C6_monthly.tif ")
|
|
| 69 | 32 |
|
| 70 |
## nobs |
|
| 71 |
#system("/usr/local/gdal-1.10.0/bin/gdalbuildvrt -allow_projection_difference -sd 4 data/MOD35C6_monthlypmiss.vrt `find /home/adamw/acrobates/adamw/projects/interp/data/modis/mod35/summary/ -name '*[0-9].nc'` ")
|
|
| 72 |
#system("/usr/local/gdal-1.10.0/bin/gdal_translate -stats -co \"COMPRESS=LZW\" -of GTiff data/MOD35C6_monthlypmiss.vrt data/MOD35C6_monthlypmiss.tif ")
|
|
| 33 |
## landcover |
|
| 34 |
if(!file.exists("data/MCD12Q1_IGBP_2009_v5_wgs84_1km.tif")){
|
|
| 35 |
system(paste("/usr/local/gdal-1.10.0/bin/gdalwarp -tr 0.008983153 0.008983153 -r mode -ot Byte -co \"COMPRESS=LZW\"",
|
|
| 36 |
" ~/acrobatesroot/jetzlab/Data/environ/global/landcover/MODIS/MCD12Q1_IGBP_2009_v5_wgs84.tif ", |
|
| 37 |
" -t_srs \"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs\" ", |
|
| 38 |
"data/MCD12Q1_IGBP_2009_v5_wgs84_1km.tif -overwrite ",sep=""))} |
|
| 39 |
lulc=raster("data/MCD12Q1_IGBP_2009_v5_wgs84_1km.tif")
|
|
| 40 |
|
|
| 41 |
## read it in |
|
| 42 |
lulc=raster("data/MCD12Q1_IGBP_2005_v51_1km_wgs84.tif")
|
|
| 43 |
# lulc=ratify(lulc) |
|
| 44 |
data(worldgrids_pal) #load palette |
|
| 45 |
IGBP=data.frame(ID=0:16,col=worldgrids_pal$IGBP[-c(18,19)], |
|
| 46 |
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)
|
|
| 47 |
IGBP$class=rownames(IGBP);rownames(IGBP)=1:nrow(IGBP) |
|
| 48 |
levels(lulc)=list(IGBP) |
|
| 49 |
lulc=crop(lulc,mod09) |
|
| 50 |
names(lulc)="MCD12Q1" |
|
| 51 |
|
|
| 52 |
## make land mask |
|
| 53 |
if(!file.exists("data/land.tif"))
|
|
| 54 |
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)
|
|
| 55 |
land=raster("data/land.tif")
|
|
| 56 |
|
|
| 73 | 57 |
|
| 74 |
#mod35c6=raster("")
|
|
| 75 |
|
|
| 76 | 58 |
## mask cloud masks to land pixels |
| 77 | 59 |
#mod09l=mask(mod09,land) |
| 78 | 60 |
#mod35l=mask(mod35,land) |
| ... | ... | |
| 81 | 63 |
### compare MOD43 and MOD17 products |
| 82 | 64 |
|
| 83 | 65 |
## MOD17 |
| 84 |
mod17=raster("data/MOD17A3_Science_NPP_mean_00_12.tif",format="GTiff")
|
|
| 85 |
NAvalue(mod17)=65535 |
|
| 86 | 66 |
#extent(mod17)=alignExtent(mod17,mod09) |
| 67 |
if(!file.exists("data/MOD17.tif"))
|
|
| 68 |
system("align.sh ~/acrobates/adamw/projects/interp/data/modis/MOD17/MOD17A3_Science_NPP_mean_00_12.tif data/MOD09_2009.tif data/MOD17.tif")
|
|
| 69 |
mod17=raster("data/MOD17.tif",format="GTiff")
|
|
| 70 |
NAvalue(mod17)=65535 |
|
| 87 | 71 |
mod17=crop(mod17,mod09) |
| 88 | 72 |
names(mod17)="MOD17" |
| 89 | 73 |
|
| 90 |
mod17qc=raster("data/MOD17A3_Science_NPP_Qc_mean_00_12.tif",format="GTiff")
|
|
| 74 |
if(!file.exists("data/MOD17qc.tif"))
|
|
| 75 |
system("align.sh ~/acrobates/adamw/projects/interp/data/modis/MOD17/MOD17A3_Science_NPP_Qc_mean_00_12.tif data/MOD09_2009.tif data/MOD17qc.tif")
|
|
| 76 |
mod17qc=raster("data/MOD17qc.tif",format="GTiff")
|
|
| 91 | 77 |
NAvalue(mod17qc)=255 |
| 92 |
#extent(mod17qc)=alignExtent(mod17qc,mod09) |
|
| 93 | 78 |
mod17qc=crop(mod17qc,mod09) |
| 94 | 79 |
names(mod17qc)="MOD17CF" |
| 95 | 80 |
|
| 96 | 81 |
## MOD11 via earth engine |
| 97 |
mod11=raster("data/MOD11_LST_2009.tif",format="GTiff")
|
|
| 82 |
if(!file.exists("data/MOD11_2009.tif"))
|
|
| 83 |
system("align.sh ~/acrobates/adamw/projects/interp/data/modis/mod11/2009/MOD11_LST_2009.tif data/MOD09_2009.tif data/MOD11_2009.tif")
|
|
| 84 |
mod11=raster("data/MOD11_2009.tif",format="GTiff")
|
|
| 98 | 85 |
names(mod11)="MOD11" |
| 99 | 86 |
NAvalue(mod11)=0 |
| 100 |
mod11qc=raster("data/MOD11_Pmiss_2009.tif",format="GTiff")
|
|
| 87 |
if(!file.exists("data/MOD11qc_2009.tif"))
|
|
| 88 |
system("align.sh ~/acrobates/adamw/projects/interp/data/modis/mod11/2009/MOD11_Pmiss_2009.tif data/MOD09_2009.tif data/MOD11qc_2009.tif")
|
|
| 89 |
mod11qc=raster("data/MOD11qc_2009.tif",format="GTiff")
|
|
| 90 |
mod11qc=crop(mod11qc,mod09) |
|
| 101 | 91 |
names(mod11qc)="MOD11CF" |
| 102 | 92 |
|
| 103 |
## MOD43 via earth engine |
|
| 104 |
#mod43=raster("data/mod43_2009.tif",format="GTiff")
|
|
| 105 |
#mod43qc=raster("data/mod43_2009.tif",format="GTiff")
|
|
| 106 |
|
|
| 107 | 93 |
|
| 108 | 94 |
### Create some summary objects for plotting |
| 109 | 95 |
#difm=v6m-v5m |
| ... | ... | |
| 111 | 97 |
#names(v5v6compare)=c("Collection 5","Collection 6","Difference (C6-C5)")
|
| 112 | 98 |
|
| 113 | 99 |
### Processing path |
| 114 |
pp=raster("data/MOD35_ProcessPath.tif")
|
|
| 115 |
#rclmat=matrix( c(0, 0, 0, 1, 1, 11, 2, 2, 16, 3, 3, 1), ncol=3, byrow=TRUE) |
|
| 116 |
#reclassify(pp,rclmat,file="data/MOD35_ProcessPath_reclass.tif",options=c("COMPRESS=LZW","ZLEVEL=9","PREDICTOR=2"),datatype="INT1U",overwrite=T)
|
|
| 117 |
#pp=raster("data/MOD35_ProcessPath_reclass.tif")
|
|
| 100 |
if(!file.exists("data/MOD35pp.tif"))
|
|
| 101 |
system("align.sh data/MOD35_ProcessPath.tif data/MOD09_2009.tif data/MOD35pp.tif")
|
|
| 102 |
pp=raster("data/MOD35pp.tif")
|
|
| 118 | 103 |
NAvalue(pp)=255 |
| 119 | 104 |
names(pp)="MOD35pp" |
| 120 |
#extent(pp)=alignExtent(pp,mod09) |
|
| 121 |
#pp=crop(pp,mod09) |
|
| 122 |
|
|
| 123 |
## Summary plot of mod17 and mod43 |
|
| 124 |
#modprod=stack(mod35c5,mod09,pp,lulc)#,mod43,mod43qc) |
|
| 125 |
#names(modprod)=c("MOD17","MOD17qc")#,"MOD43","MOD43qc")
|
|
| 126 |
|
|
| 105 |
pp=crop(pp,mod09) |
|
| 127 | 106 |
|
| 128 | 107 |
## comparison of % cloudy days |
| 129 | 108 |
dif_c5_09=mod35c5-mod09 |
| 130 | 109 |
dif_c6_09=mod35c6-mod09 |
| 131 | 110 |
dif_c5_c6=mod35c5-mod35c6 |
| 132 | 111 |
|
| 133 |
|
|
| 134 | 112 |
#hist(dif,maxsamp=1000000) |
| 135 |
|
|
| 136 | 113 |
## draw lulc-stratified random sample of mod35-mod09 differences |
| 137 | 114 |
#samp=sampleStratified(lulc, 1000, exp=10) |
| 138 | 115 |
#save(samp,file="LULC_StratifiedSample_10000.Rdata") |
| ... | ... | |
| 153 | 130 |
|
| 154 | 131 |
|
| 155 | 132 |
### Transects |
| 156 |
#r1=Lines(list( |
|
| 157 |
# Line(matrix(c( |
|
| 158 |
# -61.183,1.165, |
|
| 159 |
# -60.881,0.825 |
|
| 160 |
# -60.106, 3.809, |
|
| 161 |
# -60.847,3.240 |
|
| 162 |
# ),ncol=2,byrow=T))),"Venezuela") |
|
| 163 | 133 |
r1=Lines(list( |
| 164 | 134 |
Line(matrix(c( |
| 165 | 135 |
-61.688,4.098, |
| ... | ... | |
| 337 | 307 |
lim=c(-15,100))), |
| 338 | 308 |
xlab="Distance Along Transect (km)", ylab="% Missing Data / % of Maximum Value", |
| 339 | 309 |
legend=list( |
| 340 |
#inside=list(fun=draw.key(list( |
|
| 341 |
# lines=list(type="b",pch=16,cex=.5,lty=c(1,1,1,5,1,5),col=c("red","blue","darkgreen","darkgreen","orange","orange")),
|
|
| 342 |
# text=list(c("MOD09 % Cloudy","MOD35 % Cloudy","MOD17 % Missing","MOD17 (scaled)","MOD11 % Missing","MOD11 (scaled)"),
|
|
| 343 |
# lwd=1,col=c("red","blue","darkgreen","darkgreen","orange","orange")))),corner=c(1,.34)) ,
|
|
| 344 |
right=list(fun=draw.key(list(columns=1,#title="MCD12Q1", |
|
| 345 |
lines=list(type=c("b","b","b","b","l","b","l"),pch=16,cex=.5,lty=c(1,1,1,1,5,1,5),col=c("red","blue","black","darkgreen","darkgreen","orange","orange",rep(NA,12))),
|
|
| 346 |
rectangles=list(border=NA,col=c(rep(NA,8),"tan","darkgreen",NA,IGBP$col[sort(unique(transd$value[transd$variable=="MCD12Q1"]+1))]) ), |
|
| 347 |
text=list(c("MOD09 % Cloudy","C5 MOD35 % Cloudy","C6 MOD35 % Cloudy","MOD17 % Missing","MOD17 (scaled)","MOD11 % Missing","MOD11 (scaled)",
|
|
| 348 |
"\n \n \n C5 MOD35 Processing Path","Desert","Land","\n \n \n MCD12Q1 IGBP Land Cover",IGBP$class[sort(unique(transd$value[transd$variable=="MCD12Q1"]+1))])))))), |
|
| 349 |
strip = strip.custom(par.strip.text=list(cex=.75))) |
|
| 350 |
# |
|
| 310 |
bottom=list(fun=draw.key(list( rep=FALSE,columns=1,title=" ", |
|
| 311 |
## text=list(c("MOD09 % Cloudy","C5 MOD35 % Cloudy","C6 MOD35 % Cloudy","MOD17 % Missing","MOD17 (scaled)","MOD11 % Missing","MOD11 (scaled)")),
|
|
| 312 |
lines=list(type=c("b","b","b","b","l","b","l"),pch=16,cex=.5,lty=c(1,1,1,1,5,1,5),col=c("red","blue","black","darkgreen","darkgreen","orange","orange")),
|
|
| 313 |
text=list(c("MOD09 % Cloudy","C5 MOD35 % Cloudy","C6 MOD35 % Cloudy","MOD17 % Missing","MOD17 (scaled)","MOD11 % Missing","MOD11 (scaled)")),
|
|
| 314 |
rectangles=list(border=NA,col=c(NA,"tan","darkgreen")), |
|
| 315 |
text=list(c("C5 MOD35 Processing Path","Desert","Land")),
|
|
| 316 |
rectangles=list(border=NA,col=c(NA,IGBP$col[sort(unique(transd$value[transd$variable=="MCD12Q1"]+1))])), |
|
| 317 |
text=list(c("MCD12Q1 IGBP Land Cover",IGBP$class[sort(unique(transd$value[transd$variable=="MCD12Q1"]+1))])))))),
|
|
| 318 |
strip = strip.custom(par.strip.text=list(cex=.75))) |
|
| 351 | 319 |
print(p4) |
| 352 | 320 |
|
| 321 |
trdw=cast(trd,trans+x+y~variable,value="value") |
|
| 322 |
transdw=cast(transd,transect+dist~variable,value="value") |
|
| 323 |
|
|
| 324 |
xyplot(MOD11CF~C5MOD35CF|transect,groups=MCD12Q1,data=transdw,pch=16,cex=1,scales=list(relation="free")) |
|
| 325 |
xyplot(MOD17~C5MOD35CF|trans,groups=MCD12Q1,data=trdw,pch=16,cex=1,scales=list(relation="free")) |
|
| 326 |
|
|
| 353 | 327 |
#p5=levelplot(value~x*y|variable,data=rondonia,asp=1,scales=list(draw=F,rot=0,relation="free"),colorkey=T)#, |
| 354 | 328 |
#print(p5) |
| 355 | 329 |
|
Also available in: Unified diff
Updated profile plots of transect figure