| 1 |
99b3508d
|
Adam M. Wilson
|
############################
|
| 2 |
|
|
#### Extract MOD35 C6 processing path
|
| 3 |
|
|
|
| 4 |
d39ab57e
|
Adam M. Wilson
|
setwd("~/acrobates/adamw/projects/interp/data/modis/mod35/processpath")
|
| 5 |
99b3508d
|
Adam M. Wilson
|
library(multicore)
|
| 6 |
|
|
library(raster)
|
| 7 |
|
|
library(spgrass6)
|
| 8 |
ddd9a810
|
Adam M. Wilson
|
library(rgeos)
|
| 9 |
99b3508d
|
Adam M. Wilson
|
##download 3 days of modis swath data:
|
| 10 |
|
|
|
| 11 |
|
|
url="ftp://ladsweb.nascom.nasa.gov/allData/51/MOD35_L2/2012/"
|
| 12 |
|
|
dir.create("swath")
|
| 13 |
|
|
|
| 14 |
d39ab57e
|
Adam M. Wilson
|
getdata=F
|
| 15 |
|
|
if(getdata)
|
| 16 |
|
|
system(paste("wget -S --recursive --no-parent --no-directories -N -P swath --accept \"hdf\" --accept \"002|003|004\" ",url))
|
| 17 |
99b3508d
|
Adam M. Wilson
|
|
| 18 |
|
|
|
| 19 |
|
|
### make global raster that aligns with MODLAND tiles
|
| 20 |
|
|
## get MODLAND tile to serve as base
|
| 21 |
|
|
#system("wget http://e4ftl01.cr.usgs.gov/MOLT/MOD13A3.005/2000.02.01/MOD13A3.A2000032.h00v08.005.2006271174446.hdf")
|
| 22 |
|
|
#t=raster(paste("HDF4_EOS:EOS_GRID:\"",getwd(),"/MOD13A3.A2000032.h00v08.005.2006271174446.hdf\":MOD_Grid_monthly_1km_VI:1 km monthly NDVI",sep=""))
|
| 23 |
d39ab57e
|
Adam M. Wilson
|
t=raster(paste("../../MOD17/MOD17A3_Science_NPP_mean_00_12.tif",sep=""))
|
| 24 |
ddd9a810
|
Adam M. Wilson
|
projection(t)
|
| 25 |
99b3508d
|
Adam M. Wilson
|
|
| 26 |
|
|
## make global extent
|
| 27 |
|
|
pmodis="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs"
|
| 28 |
|
|
|
| 29 |
|
|
glb=t
|
| 30 |
ddd9a810
|
Adam M. Wilson
|
#values(glb)=NA
|
| 31 |
99b3508d
|
Adam M. Wilson
|
glb=extend(glb,extent(-180,180,-90,90))
|
| 32 |
|
|
|
| 33 |
|
|
#glb=raster(glb,crs="+proj=longlat +datum=WGS84",nrows=42500,ncols=85000)
|
| 34 |
|
|
#extent(glb)=alignExtent(projectRaster(glb,crs=projection(t),over=T),t)
|
| 35 |
|
|
#res(glb)=c(926.6254,926.6264)
|
| 36 |
|
|
#projection(glb)=pmodis
|
| 37 |
|
|
|
| 38 |
|
|
## confirm extent
|
| 39 |
|
|
#projectExtent(glb,crs="+proj=longlat +datum=WGS84")
|
| 40 |
|
|
|
| 41 |
|
|
|
| 42 |
|
|
#### Grid and mosaic the swath data
|
| 43 |
|
|
|
| 44 |
5f212fe8
|
Adam M. Wilson
|
stitch="sudo MRTDATADIR=\"/usr/local/heg/2.12/data\" PGSHOME=/usr/local/heg/2.12/TOOLKIT_MTD PWD=/home/adamw /usr/local/heg/2.12/bin/swtif"
|
| 45 |
d39ab57e
|
Adam M. Wilson
|
#stitch="/usr/local/heg/2.12/bin/swtif"
|
| 46 |
99b3508d
|
Adam M. Wilson
|
|
| 47 |
9b5125b9
|
Adam M. Wilson
|
#stitch="sudo MRTDATADIR=\"/usr/local/heg/2.11/data\" PGSHOME=/usr/local/heg/2.11/TOOLKIT_MTD PWD=/home/adamw /usr/local/heg/2.11/bin/swtif"
|
| 48 |
d39ab57e
|
Adam M. Wilson
|
files=paste(getwd(),"/",list.files("swath",pattern="hdf$",full=T),sep="")
|
| 49 |
99b3508d
|
Adam M. Wilson
|
|
| 50 |
|
|
## vars to process
|
| 51 |
|
|
vars=as.data.frame(matrix(c(
|
| 52 |
5f212fe8
|
Adam M. Wilson
|
"Cloud_Mask", "CM", "NN", 1,
|
| 53 |
|
|
# "Sensor_Azimuth", "ZA", "CUBIC", 1,
|
| 54 |
|
|
"Sensor_Zenith", "SZ", "CUBIC", 1),
|
| 55 |
|
|
byrow=T,ncol=4,dimnames=list(1:2,c("variable","varid","method","band"))),stringsAsFactors=F)
|
| 56 |
ddd9a810
|
Adam M. Wilson
|
|
| 57 |
|
|
## global bounding box
|
| 58 |
5f212fe8
|
Adam M. Wilson
|
gbb=cbind(c(-180,-180,180,180,-180),c(-90,90,90,-90,-90))
|
| 59 |
ddd9a810
|
Adam M. Wilson
|
gpp = SpatialPolygons(list(Polygons(list(Polygon(gbb)),1)))
|
| 60 |
|
|
proj4string(gpp)=projection(glb)
|
| 61 |
|
|
|
| 62 |
d39ab57e
|
Adam M. Wilson
|
outdir="~/acrobates/adamw/projects/interp/data/modis/mod35/processpath/gridded/"
|
| 63 |
ddd9a810
|
Adam M. Wilson
|
|
| 64 |
5f212fe8
|
Adam M. Wilson
|
swtif<-function(file,var){
|
| 65 |
|
|
outfile=paste(tempdir(),"/",var$varid,"_",basename(file),sep="") #gridded path
|
| 66 |
ddd9a810
|
Adam M. Wilson
|
## First write the parameter file (careful, heg is very finicky!)
|
| 67 |
5f212fe8
|
Adam M. Wilson
|
hdr=paste("NUM_RUNS = 1")
|
| 68 |
ddd9a810
|
Adam M. Wilson
|
grp=paste("
|
| 69 |
99b3508d
|
Adam M. Wilson
|
BEGIN
|
| 70 |
|
|
INPUT_FILENAME=",file,"
|
| 71 |
|
|
OBJECT_NAME=mod35
|
| 72 |
5f212fe8
|
Adam M. Wilson
|
FIELD_NAME=",var$variable,"|
|
| 73 |
|
|
BAND_NUMBER = ",var$band,"
|
| 74 |
99b3508d
|
Adam M. Wilson
|
OUTPUT_PIXEL_SIZE_X=0.008333333
|
| 75 |
|
|
OUTPUT_PIXEL_SIZE_Y=0.008333333
|
| 76 |
5f212fe8
|
Adam M. Wilson
|
SPATIAL_SUBSET_UL_CORNER = ( ",bbox(gpp)[2,2]," ",bbox(gpp)[1,1]," )
|
| 77 |
|
|
SPATIAL_SUBSET_LR_CORNER = ( ",bbox(gpp)[2,1]," ",bbox(gpp)[1,2]," )
|
| 78 |
|
|
RESAMPLING_TYPE =",var$method,"
|
| 79 |
99b3508d
|
Adam M. Wilson
|
OUTPUT_PROJECTION_TYPE = GEO
|
| 80 |
|
|
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 )
|
| 81 |
|
|
# 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 )
|
| 82 |
|
|
# projection parameters from http://landweb.nascom.nasa.gov/cgi-bin/QA_WWW/newPage.cgi?fileName=sn_gctp
|
| 83 |
|
|
ELLIPSOID_CODE = WGS84
|
| 84 |
|
|
OUTPUT_TYPE = HDFEOS
|
| 85 |
5f212fe8
|
Adam M. Wilson
|
OUTPUT_FILENAME= ",outfile,"
|
| 86 |
99b3508d
|
Adam M. Wilson
|
END
|
| 87 |
|
|
|
| 88 |
|
|
",sep="")
|
| 89 |
ddd9a810
|
Adam M. Wilson
|
## write it to a file
|
| 90 |
|
|
cat( c(hdr,grp) , file=paste(tempdir(),"/",basename(file),"_MODparms.txt",sep=""))
|
| 91 |
|
|
## now run the swath2grid tool
|
| 92 |
|
|
## write the gridded file
|
| 93 |
|
|
print(paste("Starting",file))
|
| 94 |
d39ab57e
|
Adam M. Wilson
|
system(paste("sudo ",stitch," -p ",tempdir(),"/",basename(file),"_MODparms.txt -d -log /dev/null ",sep=""),intern=F)#,ignore.stderr=F)
|
| 95 |
5f212fe8
|
Adam M. Wilson
|
print(paste("Finished processing variable",var$variable,"from ",basename(file),"to",outfile))
|
| 96 |
|
|
}
|
| 97 |
ddd9a810
|
Adam M. Wilson
|
|
| 98 |
5f212fe8
|
Adam M. Wilson
|
getpath<- function(file){
|
| 99 |
|
|
setwd(tempdir())
|
| 100 |
|
|
bfile=sub(".hdf","",basename(file))
|
| 101 |
|
|
tempfile2_path=paste(tempdir(),"/",bfile,".tif",sep="") #gridded/masked/processed path
|
| 102 |
|
|
outfile=paste(outdir,"/",bfile,".tif",sep="") #final file
|
| 103 |
d39ab57e
|
Adam M. Wilson
|
if(file.exists(outfile)) return(paste(file," already finished"))
|
| 104 |
5f212fe8
|
Adam M. Wilson
|
ppc=gpp
|
| 105 |
|
|
#######
|
| 106 |
|
|
## run swtif for each band
|
| 107 |
|
|
lapply(1:nrow(vars),function(i) swtif(file,vars[i,]))
|
| 108 |
ddd9a810
|
Adam M. Wilson
|
####### import to R for processing
|
| 109 |
5f212fe8
|
Adam M. Wilson
|
|
| 110 |
|
|
if(!file.exists(paste(tempdir(),"/CM_",basename(file),sep=""))) {
|
| 111 |
ddd9a810
|
Adam M. Wilson
|
file.remove(list.files(tempdir(),pattern=bfile,full=T))
|
| 112 |
|
|
return(c(file,0))
|
| 113 |
|
|
}
|
| 114 |
|
|
## convert to land path
|
| 115 |
5f212fe8
|
Adam M. Wilson
|
d=raster(paste(tempdir(),"/CM_",basename(file),sep=""))
|
| 116 |
|
|
sz=raster(paste(tempdir(),"/SZ_",basename(file),sep=""))
|
| 117 |
|
|
NAvalue(sz)=-9999
|
| 118 |
d39ab57e
|
Adam M. Wilson
|
getlc=function(x,y) {ifelse(y<0|y>6000,NA,((x%/%2^6) %% 2^2))}
|
| 119 |
ddd9a810
|
Adam M. Wilson
|
path= overlay(d,sz,fun=getlc,filename=tempfile2_path,options=c("COMPRESS=LZW", "LEVEL=9","PREDICTOR=2"),datatype="INT1U",overwrite=T)
|
| 120 |
|
|
### warp them to align all pixels
|
| 121 |
|
|
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=""))
|
| 122 |
|
|
## delete temporary files
|
| 123 |
|
|
file.remove(list.files(tempdir(),pattern=bfile,full=T))
|
| 124 |
|
|
return(c(file,1))
|
| 125 |
|
|
}
|
| 126 |
|
|
|
| 127 |
|
|
|
| 128 |
5f212fe8
|
Adam M. Wilson
|
### run it
|
| 129 |
905db3dc
|
Adam M. Wilson
|
mclapply(files,getpath,mc.cores=10)
|
| 130 |
99b3508d
|
Adam M. Wilson
|
|
| 131 |
|
|
## check gdal can read all of them
|
| 132 |
5f212fe8
|
Adam M. Wilson
|
gfiles=list.files(outdir,pattern="tif$",full=T)
|
| 133 |
99b3508d
|
Adam M. Wilson
|
length(gfiles)
|
| 134 |
|
|
|
| 135 |
|
|
check=do.call(rbind,mclapply(gfiles,function(file){
|
| 136 |
|
|
gd=system(paste("gdalinfo ",file,sep=""),intern=T)
|
| 137 |
|
|
if(any(grepl("Corner",gd))) return(1)
|
| 138 |
|
|
else return(0)
|
| 139 |
|
|
}))
|
| 140 |
|
|
|
| 141 |
|
|
table(check)
|
| 142 |
|
|
|
| 143 |
|
|
file.remove(gfiles[check==0])
|
| 144 |
|
|
|
| 145 |
ddd9a810
|
Adam M. Wilson
|
## use new gdal
|
| 146 |
d39ab57e
|
Adam M. Wilson
|
#system(paste("nohup /usr/local/gdal-1.10.0/bin/gdalwarp -wm 900 -overwrite -co COMPRESS=LZW -co PREDICTOR=2 -multi ",outdir,"/*.tif MOD35_path_gdalwarp.tif &",sep=""))
|
| 147 |
|
|
#system(paste("nohup /usr/local/gdal-1.10.0/bin/gdalwarp -wm 900 -overwrite -co COMPRESS=LZW -co PREDICTOR=2 -multi -r mode ",outdir,"/*.tif MOD35_path_gdalwarp_mode.tif &",sep=""))
|
| 148 |
|
|
#system(paste(" /usr/local/gdal-1.10.0/bin/gdalwarp -wm 900 -overwrite -co COMPRESS=LZW -co PREDICTOR=2 -multi ",outdir,"/MOD35_L2.A2012001*.tif MOD35_path_gdalwarp.tif",sep=""))
|
| 149 |
99b3508d
|
Adam M. Wilson
|
|
| 150 |
|
|
|
| 151 |
be7b1081
|
Adam M. Wilson
|
### Merge them into a geotiff
|
| 152 |
d39ab57e
|
Adam M. Wilson
|
#system(paste("/usr/local/src/gdal-1.10.0/swig/python/scripts/gdal_merge.py -v -init 255 -n 255 -a_nodata 255 -o MOD35_ProcessPath_gdalmerge2.tif -co \"ZLEVEL=9\" -co \"COMPRESS=LZW\" -co \"PREDICTOR=2\" `ls -d -1 ",outdir,"/*.tif --sort=size | head -n 20 ` ",sep=""))
|
| 153 |
|
|
system(paste("/usr/local/src/gdal-1.10.0/swig/python/scripts/gdal_merge.py -v -o MOD35_ProcessPath_gdalmerge2.tif -co \"ZLEVEL=9\" -co \"COMPRESS=LZW\" -co \"PREDICTOR=2\" `ls -d -1 ",outdir,"/*.tif --sort=size ` &",sep=""))
|
| 154 |
be7b1081
|
Adam M. Wilson
|
|
| 155 |
d39ab57e
|
Adam M. Wilson
|
|
| 156 |
99b3508d
|
Adam M. Wilson
|
## try with pktools
|
| 157 |
|
|
## global
|
| 158 |
d39ab57e
|
Adam M. Wilson
|
#system(paste("pkmosaic -co COMPRESS=LZW -co PREDICTOR=2 ",paste("-i",list.files("gridded",full=T,pattern="tif$")[1:10],collapse=" ")," -o MOD35_path_pkmosaic_mode.tif -m 6 -v -t 255 -t 0 &"))
|
| 159 |
99b3508d
|
Adam M. Wilson
|
#bb="-ulx -180 -uly 90 -lrx 180 -lry -90"
|
| 160 |
|
|
#bb="-ulx -180 -uly 90 -lrx 170 -lry 80"
|
| 161 |
d39ab57e
|
Adam M. Wilson
|
#bb="-ulx -72 -uly 11 -lrx -59 -lry -1"
|
| 162 |
99b3508d
|
Adam M. Wilson
|
|
| 163 |
|
|
|
| 164 |
|
|
#expand.grid(x=seq(-180,170,by=10),y=seq(-90,80))
|
| 165 |
d39ab57e
|
Adam M. Wilson
|
#gf2= grep("2012009[.]03",gfiles,value=T)
|
| 166 |
|
|
#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"))
|
| 167 |
99b3508d
|
Adam M. Wilson
|
|
| 168 |
|
|
# bounding box?
|
| 169 |
|
|
|
| 170 |
|
|
###########
|
| 171 |
|
|
### Use GRASS to import all the tifs and calculat the mode
|
| 172 |
|
|
## make temporary working directory
|
| 173 |
|
|
tf=paste(tempdir(),"/grass", Sys.getpid(),"/", sep="") #temporar
|
| 174 |
|
|
if(!file.exists(tf)) dir.create(tf)
|
| 175 |
|
|
|
| 176 |
|
|
## set up temporary grass instance for this PID
|
| 177 |
|
|
gisBase="/usr/lib/grass64"
|
| 178 |
|
|
print(paste("Set up temporary grass session in",tf))
|
| 179 |
|
|
initGRASS(gisBase=gisBase,gisDbase=tf,SG=as(glb,"SpatialGridDataFrame"),override=T,location="mod35",mapset="PERMANENT",home=tf,pid=Sys.getpid())
|
| 180 |
|
|
system(paste("g.proj -c proj4=\"",projection(glb),"\"",sep=""),ignore.stdout=T,ignore.stderr=T)
|
| 181 |
|
|
|
| 182 |
|
|
## read in NPP grid to serve as grid
|
| 183 |
|
|
execGRASS("r.in.gdal",input=t@file@name,output="grid")
|
| 184 |
|
|
system("g.region rast=grid n=90 s=-90 save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)
|
| 185 |
|
|
|
| 186 |
|
|
## get files already imported - only important if more tifs were written after an initial import
|
| 187 |
|
|
imported=basename(gfiles)%in%system("g.mlist type=rast pattern=MOD*",intern=T)
|
| 188 |
|
|
table(imported)
|
| 189 |
|
|
|
| 190 |
|
|
## read in all tifs
|
| 191 |
d39ab57e
|
Adam M. Wilson
|
for(f in gfiles[!imported]) {
|
| 192 |
|
|
print(f)
|
| 193 |
|
|
execGRASS("r.in.gdal",input=f,output=basename(f),flags="o")
|
| 194 |
|
|
}
|
| 195 |
|
|
|
| 196 |
|
|
## calculate mode - can't have more than 1000 open files
|
| 197 |
|
|
execGRASS("r.series",input=paste(system("g.mlist type=rast pattern=MOD*",intern=T)[1:1000],sep="",collapse=","),output="path1",method="mode",range=c(1,5),flags=c("verbose","overwrite"))
|
| 198 |
|
|
execGRASS("r.series",input=paste(system("g.mlist type=rast pattern=MOD*",intern=T)[1001:2000],sep="",collapse=","),output="path2",method="mode",range=c(1,5),flags=c("verbose","overwrite"))
|
| 199 |
|
|
execGRASS("r.series",input=paste(system("g.mlist type=rast pattern=MOD*",intern=T)[2001:3000],sep="",collapse=","),output="path3",method="mode",range=c(1,5),flags=c("verbose","overwrite"))
|
| 200 |
|
|
execGRASS("r.series",input=paste(system("g.mlist type=rast pattern=MOD*",intern=T)[3001:4000],sep="",collapse=","),output="path4",method="mode",range=c(1,5),flags=c("verbose","overwrite"))
|
| 201 |
|
|
execGRASS("r.series",input=paste(system("g.mlist type=rast pattern=MOD*",intern=T)[4001:5000],sep="",collapse=","),output="path5",method="mode",range=c(1,5),flags=c("verbose","overwrite"))
|
| 202 |
|
|
|
| 203 |
|
|
## Get mode of modes
|
| 204 |
|
|
execGRASS("r.series",input=paste(system("g.mlist type=rast pattern=path*",intern=T),sep="",collapse=","),output="MOD35_path",method="mode",range=c(1,5),flags=c("verbose","overwrite"))
|
| 205 |
99b3508d
|
Adam M. Wilson
|
|
| 206 |
|
|
## add colors
|
| 207 |
|
|
execGRASS("r.colors",map="MOD35_path",rules="MOD35_path_grasscolors.txt")
|
| 208 |
|
|
## write to disk
|
| 209 |
|
|
execGRASS("r.out.gdal",input="MOD35_path",output=paste(getwd(),"/MOD35_ProcessPath_C5.tif",sep=""),type="Byte",createopt="COMPRESS=LZW,LEVEL=9,PREDICTOR=2")
|
| 210 |
|
|
|
| 211 |
|
|
### delete the temporary files
|
| 212 |
|
|
unlink_.gislock()
|
| 213 |
|
|
system(paste("rm -frR ",tf,sep=""))
|
| 214 |
|
|
|
| 215 |
|
|
#########################
|
| 216 |
|
|
|
| 217 |
|
|
|
| 218 |
|
|
cols=c("blue","lightblue","tan","green")
|
| 219 |
|
|
|
| 220 |
|
|
|
| 221 |
|
|
|
| 222 |
|
|
## connect to raster to extract land-cover bit
|
| 223 |
|
|
library(raster)
|
| 224 |
|
|
|
| 225 |
|
|
d=raster("CM.tif")
|
| 226 |
|
|
getlc=function(x) {(x/2^6) %% 2^2}
|
| 227 |
|
|
|
| 228 |
|
|
calc(d,fun=getlc,filename="CM_LC.tif")
|
