Project

General

Profile

« Previous | Next » 

Revision 58d05081

Added by Adam Wilson about 11 years ago

Adding ee.MOD09.py script to process and download data from Google Earth Engine

View differences:

climate/procedures/MOD35C5_Evaluation.r
11 11 
library(rgeos)
12 12 
library(splancs)
13 13 

  		




  
  14
  
    

  		




  
  15
  
    
## add tags for distribution


  		




  
  16
  
    
## MOD35


  		




  
  17
  
    
tags=c("TIFFTAG_IMAGEDESCRIPTION='Collection 5 MOD35 Cloud Frequency for 2009 extracted from MOD09GA state_1km bits 0-1. The MOD35 bits encode four categories (with associated confidence that the pixel is actually clear): confidently clear (confidence > 0.99), probably clear (0.99 >= confidence > 0.95), probably cloudy (0.95 >= confidence > 0.66), and confidently cloudy (confidence <= 0.66).  Following the advice of the MODIS science team (Frey, 2010), we binned confidently clear and probably clear together as clear and the other two classes as cloudy.  The daily cloud mask time series were summarized to mean cloud frequency (CF) by calculating the proportion of cloudy days during 2009'",


  		




  
  18
  
    
  "TIFFTAG_DOCUMENTNAME='Collection 5 MOD35 Cloud Frequency'",


  		




  
  19
  
    
  "TIFFTAG_DATETIME='20090101'",


  		




  
  20
  
    
  "TIFFTAG_ARTIST='Adam M. Wilson (adam.wilson@yale.edu)'")


  		




  
  21
  
    
system(paste("/usr/local/src/gdal-1.10.0/swig/python/scripts/gdal_edit.py data/MOD35_2009.tif ",paste("-mo ",tags,sep="",collapse=" "),sep=""))


  		




  
  22
  
    
## MOD09


  		




  
  23
  
    
tags=c("TIFFTAG_IMAGEDESCRIPTION='Collection 5 MOD09 Cloud Frequency for 2009 extracted from MOD09GA \'PGE11\' internal cloud mask algorithm (embedded in MOD09GA \'state_1km\' bit 10. The daily cloud mask time series were summarized to mean cloud frequency (CF) by calculating the proportion of cloudy days during 2009'",


  		




  
  24
  
    
  "TIFFTAG_DOCUMENTNAME='Collection 5 MOD09 Cloud Frequency'",


  		




  
  25
  
    
  "TIFFTAG_DATETIME='20090101'",


  		




  
  26
  
    
  "TIFFTAG_ARTIST='Adam M. Wilson (adam.wilson@yale.edu)'")


  		




  
  27
  
    
system(paste("/usr/local/src/gdal-1.10.0/swig/python/scripts/gdal_edit.py data/MOD09_2009.tif ",paste("-mo ",tags,sep="",collapse=" "),sep=""))


  		




  
  28
  
    

  		




  14
  29
  
    
## get % cloudy


  		




  15
  30
  
    
mod09=raster("data/MOD09_2009.tif")


  		




  16
  31
  
    
names(mod09)="C5MOD09CF"


  		




  ......




  20
  35
  
    
names(mod35c5)="C5MOD35CF"


  		




  21
  36
  
    
NAvalue(mod35c5)=0


  		




  22
  37
  
    

  		




  
  38
  
    

  		




  23
  39
  
    
## mod35C6 annual


  		




  24
  40
  
    
if(!file.exists("data/MOD35C6_2009.tif")){


  		




  25
  41
  
    
  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 `find /home/adamw/acrobates/adamw/projects/interp/data/modis/mod35/summary/ -name '*h[1-9]*_mean.nc'` ")


  		




  26
  
  
    
#  system("gdalbuildvrt data/MOD35C6.vrt `find /home/adamw/acrobates/adamw/projects/interp/data/modis/mod35/summary/ -name '*h[1-9]*_mean.nc'` ")


  		




  27
  
  
    

  		




  28
  
  
    
#  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 4 -b 1 data/MOD35C6_CFday_pmiss.vrt `find /home/adamw/acrobates/adamw/projects/interp/data/modis/mod35/summary/ -name '*h[1]*.nc'` ")


  		




  29
  
  
    
#  system("gdalwarp data/MOD35C6_CFday_pmiss.vrt data/MOD35C6_CFday_pmiss.tif -r bilinear")


  		




  30
  
  
    

  		




  31
  42
  
    
  system("align.sh data/MOD35C6.vrt data/MOD09_2009.tif data/MOD35C6_2009.tif")


  		




  32
  
  
    
#  system("align.sh data/MOD35C6_CFday_pmiss.vrt data/MOD09_2009.tif data/MOD35C6_CFday_pmiss.tif")


  		




  33
  43
  
    
}


  		




  
  44
  
    

  		




  34
  45
  
    
mod35c6=raster("data/MOD35C6_2009.tif")


  		




  35
  46
  
    
names(mod35c6)="C6MOD35CF"


  		




  36
  47
  
    
NAvalue(mod35c6)=255


  		




  ......




  44
  55
  
    
               "data/MCD12Q1_IGBP_2009_051_wgs84_1km.tif -overwrite ",sep=""))}


  		




  45
  56
  
    
lulc=raster("data/MCD12Q1_IGBP_2009_051_wgs84_1km.tif")


  		




  46
  57
  
    

  		




  47
  
  
    
#  lulc=ratify(lulc)


  		




  48
  58
  
    
  data(worldgrids_pal)  #load palette


  		




  49
  59
  
    
  IGBP=data.frame(ID=0:16,col=worldgrids_pal$IGBP[-c(18,19)],


  		




  50
  60
  
    
    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)


  		




  51
  61
  
    
  IGBP$class=rownames(IGBP);rownames(IGBP)=1:nrow(IGBP)


  		




  52
  62
  
    
  levels(lulc)=list(IGBP)


  		




  53
  
  
    
#lulc=crop(lulc,mod09)


  		




  54
  63
  
    
names(lulc)="MCD12Q1"


  		




  55
  64
  
    

  		




  56
  65
  
    
## make land mask


  		




  ......




  58
  67
  
    
  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)


  		




  59
  68
  
    
land=raster("data/land.tif")


  		




  60
  69
  
    

  		




  61
  
  
    
## mask cloud masks to land pixels


  		




  62
  
  
    
#mod09l=mask(mod09,land)


  		




  63
  
  
    
#mod35l=mask(mod35,land)


  		




  64
  
  
    

  		




  65
  70
  
    
#####################################


  		




  66
  71
  
    
### compare MOD43 and MOD17 products


  		




  67
  72
  
    

  		




  68
  73
  
    
## MOD17


  		




  69
  
  
    
#extent(mod17)=alignExtent(mod17,mod09)


  		




  70
  74
  
    
if(!file.exists("data/MOD17.tif"))


  		




  71
  75
  
    
system("align.sh ~/acrobates/adamw/projects/interp/data/modis/MOD17/MOD17A3_Science_NPP_mean_00_12.tif data/MOD09_2009.tif data/MOD17.tif")


  		




  72
  76
  
    
mod17=raster("data/MOD17.tif",format="GTiff")


  		




  ......




  98
  102
  
    
names(pp)="MOD35pp"


  		




  99
  103
  
    

  		




  100
  104
  
    

  		




  101
  
  
    
#hist(dif,maxsamp=1000000)


  		




  102
  
  
    
## draw lulc-stratified random sample of mod35-mod09 differences 


  		




  103
  
  
    
#samp=sampleStratified(lulc, 1000, exp=10)


  		




  104
  
  
    
#save(samp,file="LULC_StratifiedSample_10000.Rdata")


  		




  105
  
  
    
#mean(dif[samp],na.rm=T)


  		




  106
  
  
    
#Stats(dif,function(x) c(mean=mean(x),sd=sd(x)))


  		




  107
  
  
    

  		




  108
  
  
    

  		




  109
  105
  
    
###


  		




  110
  
  
    

  		




  111
  106
  
    
n=100


  		




  112
  107
  
    
at=seq(0,100,len=n)


  		




  113
  108
  
    
cols=grey(seq(0,1,len=n))


  		




  ......




  115
  110
  
    
bgyr=colorRampPalette(c("blue","green","yellow","red"))


  		




  116
  111
  
    
cols=bgyr(n)


  		




  117
  112
  
    

  		




  118
  
  
    

  		




  119
  113
  
    
### Transects


  		




  120
  114
  
    
r1=Lines(list(


  		




  121
  115
  
    
  Line(matrix(c(


  		




  ......




  209
  203
  
    
  overlay(mod35c5,mod09,fun=function(x,y) {return(x-y)},file="data/dif_c5_09.tif",format="GTiff",options=c("COMPRESS=LZW","ZLEVEL=9"),overwrite=T)


  		




  210
  204
  
    
dif_c5_09=raster("data/dif_c5_09.tif",format="GTiff")


  		




  211
  205
  
    

  		




  212
  
  
    
#dif_c6_09=mod35c6-mod09


  		




  213
  
  
    
#dif_c5_c6=mod35c5-mod35c6


  		




  214
  206
  
    

  		




  215
  207
  
    
##################################################################################


  		




  216
  208
  
    
## Identify problematic areas with hard edges in cloud frequency


  		




  ......




  310
  302
  
    
### read them back in


  		




  311
  303
  
    
pp_bias=raster("data/pp_bias.tif")


  		




  312
  304
  
    
names(pp_bias)="Processing Path"


  		




  
  305
  
    
NAvalue(pp_bias)=255


  		




  313
  306
  
    
lulc_bias=raster("data/lulc_bias.tif")


  		




  314
  307
  
    
names(lulc_bias)="Land Use Land Cover"


  		




  
  308
  
    
NAvalue(lulc_bias)=255


  		




  
  309
  
    

  		




  
  310
  
    
## read in WWF biome data to summarize by biome


  		




  
  311
  
    
if(!file.exists("data/teow/wwf_terr_ecos.shp"){


  		




  
  312
  
    
  system("wget -O data/teow.zip http://assets.worldwildlife.org/publications/15/files/original/official_teow.zip?1349272619")


  		




  
  313
  
    
  system("unzip data/teow.zip -d data/teow/")


  		




  
  314
  
    
   biome=readOGR("data/teow/wwf_terr_ecos.shp","wwf_terr_ecos")


  		




  
  315
  
    
   biome=biome[biome$BIOME<50,]


  		




  
  316
  
    
   biome2=gUnaryUnion(biome,id=biome$BIOME)


  		




  
  317
  
    
  ## create biome.csv using names in html file   


  		




  
  318
  
    
  biomeid=read.csv("data/teow/biome.csv",stringsAsFactors=F)


  		




  
  319
  
    
  biome2=SpatialPolygonsDataFrame(biome2,data=biomeid)


  		




  
  320
  
    
  writeOGR(biome2,"data/teow","biomes",driver="ESRI Shapefile")


  		




  
  321
  
    
}


  		




  
  322
  
    
   biome=readOGR("data/teow/biomes.shp","biomes")


  		




  
  323
  
    

  		




  
  324
  
    
biome2=extract(pp_bias,biome[biome$Biome==12,],df=T,fun=function(x) data.frame(mean=mean(x,na.rm=T),sd=sd(x,na.rm=T),prop=(sum(!is.na(x))/length(x))))


  		




  315
  325
  
    

  		




  316
  
  
    
pat=c(0,0.05,1)#seq(0,0.-5,len=2) #,seq(0.05,.1,len=50))


  		




  317
  
  
    
grayr2=colorRampPalette(c("red","transparent"))#grey(c(.75,.5,.25))))


  		




  
  326
  
    
pat=c(seq(0,100,len=100),254)#seq(0,0.-5,len=2) #,seq(0.05,.1,len=50))


  		




  
  327
  
    
grayr2=colorRampPalette(c("grey","green","red"))#


  		




  
  328
  
    
#grayr2=colorRampPalette(grey(c(.75,.5,.25)))


  		




  318
  329
  
    
levelplot(stack(pp_bias,lulc_bias),col.regions=c(grayr2(2)),at=pat,


  		




  319
  330
  
    
          colorkey=F,margin=F,maxpixels=1e6)+layer(sp.lines(coast,lwd=.5))


  		




  
  331
  
    
levelplot(lulc_bias,col.regions=c(grayr2(100),"black"),at=pat,


  		




  
  332
  
    
          colorkey=T,margin=F,maxpixels=1e4)+layer(sp.lines(coast,lwd=.5))


  		




  
  333
  
    

  		




  
  334
  
    
histogram(lulc_bias)


  		




  320
  335
  
    

  		




  321
  336
  
    
cor(td1$MOD17,td1$C6MOD35,use="complete",method="spearman")


  		




  322
  337
  
    
cor(td1$MOD17[td1$edgeb==1],td1$C5MOD35[td1$edgeb==1],use="complete",method="spearman")


  		












Also available in:  Unified diff