| 1 |
99b3508d
|
Adam M. Wilson
|
############################
|
| 2 |
|
|
#### Extract MOD35 C6 processing path
|
| 3 |
|
|
|
| 4 |
|
|
setwd("~/acrobates/adamw/projects/interp/data/modis/mod35")
|
| 5 |
|
|
library(multicore)
|
| 6 |
|
|
library(raster)
|
| 7 |
|
|
library(spgrass6)
|
| 8 |
|
|
|
| 9 |
|
|
##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 |
|
|
system(paste("wget -S --recursive --no-parent --no-directories -N -P swath --accept \"hdf\" --accept \"002|003|004\" ",url))
|
| 15 |
|
|
|
| 16 |
|
|
|
| 17 |
|
|
### make global raster that aligns with MODLAND tiles
|
| 18 |
|
|
## get MODLAND tile to serve as base
|
| 19 |
|
|
#system("wget http://e4ftl01.cr.usgs.gov/MOLT/MOD13A3.005/2000.02.01/MOD13A3.A2000032.h00v08.005.2006271174446.hdf")
|
| 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=""))
|
| 22 |
|
|
|
| 23 |
|
|
## make global extent
|
| 24 |
|
|
pmodis="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs"
|
| 25 |
|
|
|
| 26 |
|
|
glb=t
|
| 27 |
|
|
values(glb)=NA
|
| 28 |
|
|
glb=extend(glb,extent(-180,180,-90,90))
|
| 29 |
|
|
|
| 30 |
|
|
#glb=raster(glb,crs="+proj=longlat +datum=WGS84",nrows=42500,ncols=85000)
|
| 31 |
|
|
#extent(glb)=alignExtent(projectRaster(glb,crs=projection(t),over=T),t)
|
| 32 |
|
|
#res(glb)=c(926.6254,926.6264)
|
| 33 |
|
|
#projection(glb)=pmodis
|
| 34 |
|
|
|
| 35 |
|
|
## confirm extent
|
| 36 |
|
|
#projectExtent(glb,crs="+proj=longlat +datum=WGS84")
|
| 37 |
|
|
|
| 38 |
|
|
|
| 39 |
|
|
#### Grid and mosaic the swath data
|
| 40 |
|
|
|
| 41 |
|
|
stitch="sudo MRTDATADIR=\"/usr/local/heg/data\" PGSHOME=/usr/local/heg/TOOLKIT_MTD PWD=/home/adamw /usr/local/heg/bin/swtif"
|
| 42 |
|
|
|
| 43 |
|
|
files=paste(getwd(),"/",list.files("swath",pattern="hdf$",full=T),sep="")
|
| 44 |
|
|
|
| 45 |
|
|
## vars to process
|
| 46 |
|
|
vars=as.data.frame(matrix(c(
|
| 47 |
|
|
"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("
|
| 66 |
|
|
BEGIN
|
| 67 |
|
|
INPUT_FILENAME=",file,"
|
| 68 |
|
|
OBJECT_NAME=mod35
|
| 69 |
|
|
FIELD_NAME=",vars$variable,"|
|
| 70 |
|
|
BAND_NUMBER = ",1:length(vars$varid),"
|
| 71 |
|
|
OUTPUT_PIXEL_SIZE_X=0.008333333
|
| 72 |
|
|
OUTPUT_PIXEL_SIZE_Y=0.008333333
|
| 73 |
|
|
SPATIAL_SUBSET_UL_CORNER = ( 90 -180 )
|
| 74 |
|
|
SPATIAL_SUBSET_LR_CORNER = ( -90 180 )
|
| 75 |
|
|
OUTPUT_OBJECT_NAME = mod35|
|
| 76 |
|
|
RESAMPLING_TYPE =NN
|
| 77 |
|
|
OUTPUT_PROJECTION_TYPE = GEO
|
| 78 |
|
|
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 )
|
| 79 |
|
|
# 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 )
|
| 80 |
|
|
# projection parameters from http://landweb.nascom.nasa.gov/cgi-bin/QA_WWW/newPage.cgi?fileName=sn_gctp
|
| 81 |
|
|
ELLIPSOID_CODE = WGS84
|
| 82 |
|
|
OUTPUT_TYPE = HDFEOS
|
| 83 |
|
|
OUTPUT_FILENAME= ",tempfile,"
|
| 84 |
|
|
END
|
| 85 |
|
|
|
| 86 |
|
|
",sep="")
|
| 87 |
|
|
|
| 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 |
|
|
})
|
| 118 |
|
|
|
| 119 |
|
|
## check gdal can read all of them
|
| 120 |
|
|
gfiles=list.files("gridded",pattern="tif$",full=T)
|
| 121 |
|
|
length(gfiles)
|
| 122 |
|
|
|
| 123 |
|
|
check=do.call(rbind,mclapply(gfiles,function(file){
|
| 124 |
|
|
gd=system(paste("gdalinfo ",file,sep=""),intern=T)
|
| 125 |
|
|
if(any(grepl("Corner",gd))) return(1)
|
| 126 |
|
|
else return(0)
|
| 127 |
|
|
}))
|
| 128 |
|
|
|
| 129 |
|
|
table(check)
|
| 130 |
|
|
|
| 131 |
|
|
file.remove(gfiles[check==0])
|
| 132 |
|
|
|
| 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"))
|
| 146 |
|
|
|
| 147 |
|
|
|
| 148 |
|
|
## try with pktools
|
| 149 |
|
|
## 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 &"))
|
| 151 |
|
|
#bb="-ulx -180 -uly 90 -lrx 180 -lry -90"
|
| 152 |
|
|
#bb="-ulx -180 -uly 90 -lrx 170 -lry 80"
|
| 153 |
|
|
bb="-ulx -72 -uly 11 -lrx -59 -lry -1"
|
| 154 |
|
|
|
| 155 |
|
|
|
| 156 |
|
|
#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 &"))
|
| 159 |
|
|
|
| 160 |
|
|
# bounding box?
|
| 161 |
|
|
|
| 162 |
|
|
###########
|
| 163 |
|
|
### Use GRASS to import all the tifs and calculat the mode
|
| 164 |
|
|
## make temporary working directory
|
| 165 |
|
|
tf=paste(tempdir(),"/grass", Sys.getpid(),"/", sep="") #temporar
|
| 166 |
|
|
if(!file.exists(tf)) dir.create(tf)
|
| 167 |
|
|
|
| 168 |
|
|
## set up temporary grass instance for this PID
|
| 169 |
|
|
gisBase="/usr/lib/grass64"
|
| 170 |
|
|
print(paste("Set up temporary grass session in",tf))
|
| 171 |
|
|
initGRASS(gisBase=gisBase,gisDbase=tf,SG=as(glb,"SpatialGridDataFrame"),override=T,location="mod35",mapset="PERMANENT",home=tf,pid=Sys.getpid())
|
| 172 |
|
|
system(paste("g.proj -c proj4=\"",projection(glb),"\"",sep=""),ignore.stdout=T,ignore.stderr=T)
|
| 173 |
|
|
|
| 174 |
|
|
## read in NPP grid to serve as grid
|
| 175 |
|
|
execGRASS("r.in.gdal",input=t@file@name,output="grid")
|
| 176 |
|
|
system("g.region rast=grid n=90 s=-90 save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)
|
| 177 |
|
|
|
| 178 |
|
|
## get files already imported - only important if more tifs were written after an initial import
|
| 179 |
|
|
imported=basename(gfiles)%in%system("g.mlist type=rast pattern=MOD*",intern=T)
|
| 180 |
|
|
table(imported)
|
| 181 |
|
|
|
| 182 |
|
|
## read in all tifs
|
| 183 |
|
|
for(f in gfiles[!imported]) execGRASS("r.in.gdal",input=f,output=basename(f),flags="o")
|
| 184 |
|
|
|
| 185 |
|
|
## calculate mode
|
| 186 |
|
|
execGRASS("r.series",input=paste(system("g.mlist type=rast pattern=MOD*",intern=T)[1:1000],sep="",collapse=","),output="MOD35_path",method="mode",range=c(1,5),flags=c("verbose","overwrite"))
|
| 187 |
|
|
## add colors
|
| 188 |
|
|
execGRASS("r.colors",map="MOD35_path",rules="MOD35_path_grasscolors.txt")
|
| 189 |
|
|
## write to disk
|
| 190 |
|
|
execGRASS("r.out.gdal",input="MOD35_path",output=paste(getwd(),"/MOD35_ProcessPath_C5.tif",sep=""),type="Byte",createopt="COMPRESS=LZW,LEVEL=9,PREDICTOR=2")
|
| 191 |
|
|
|
| 192 |
|
|
### delete the temporary files
|
| 193 |
|
|
unlink_.gislock()
|
| 194 |
|
|
system(paste("rm -frR ",tf,sep=""))
|
| 195 |
|
|
|
| 196 |
|
|
|
| 197 |
|
|
#########################
|
| 198 |
|
|
|
| 199 |
|
|
|
| 200 |
|
|
cols=c("blue","lightblue","tan","green")
|
| 201 |
|
|
|
| 202 |
|
|
|
| 203 |
|
|
### Merge them into a geotiff
|
| 204 |
|
|
system(paste("gdal_merge.py -v -n 255 -o MOD35_ProcessPath.tif -co \"COMPRESS=LZW\" -co \"PREDICTOR=2\" `ls -d -1 gridded/*.tif`",sep=""))
|
| 205 |
|
|
|
| 206 |
|
|
|
| 207 |
|
|
## connect to raster to extract land-cover bit
|
| 208 |
|
|
library(raster)
|
| 209 |
|
|
|
| 210 |
|
|
d=raster("CM.tif")
|
| 211 |
|
|
getlc=function(x) {(x/2^6) %% 2^2}
|
| 212 |
|
|
|
| 213 |
|
|
calc(d,fun=getlc,filename="CM_LC.tif")
|
