Revision 555815c9
Added by Adam Wilson about 11 years ago
| climate/procedures/MOD35_ExtractProcessPath.r | ||
|---|---|---|
| 8 | 8 |
library(rgeos) |
| 9 | 9 |
##download 3 days of modis swath data: |
| 10 | 10 |
|
| 11 |
#### Set up command for running swtif to grid and mosaic the swath data |
|
| 12 |
stitch=paste("sudo MRTDATADIR=\"/usr/local/heg/2.12/data\" PGSHOME=/usr/local/heg/2.12/TOOLKIT_MTD PWD=/home/adamw /usr/local/heg/2.12b/bin/swtif")
|
|
| 13 |
|
|
| 14 |
## Link to MOD35 Swath data |
|
| 11 | 15 |
url="ftp://ladsweb.nascom.nasa.gov/allData/51/MOD35_L2/2012/" |
| 12 | 16 |
dir.create("swath")
|
| 13 | 17 |
|
| ... | ... | |
| 15 | 19 |
if(getdata) |
| 16 | 20 |
system(paste("wget -S --recursive --no-parent --no-directories -N -P swath --accept \"hdf\" --accept \"002|003|004\" ",url))
|
| 17 | 21 |
|
| 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=""))
|
|
| 22 |
### make global raster that aligns with MOD17 NPP raster |
|
| 23 | 23 |
t=raster(paste("../../MOD17/MOD17A3_Science_NPP_mean_00_12.tif",sep=""))
|
| 24 | 24 |
projection(t) |
| 25 | 25 |
|
| 26 | 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 | 27 |
glb=t |
| 30 |
#values(glb)=NA |
|
| 31 | 28 |
glb=extend(glb,extent(-180,180,-90,90)) |
| 32 | 29 |
|
| 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 |
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 |
#stitch="/usr/local/heg/2.12/bin/swtif" |
|
| 46 |
|
|
| 47 |
#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 |
files=paste(getwd(),"/",list.files("swath",pattern="hdf$",full=T),sep="")
|
|
| 30 |
### list of swath files |
|
| 31 |
files=paste(getwd(),"/",list.files("swath",pattern="hdf$",full=T),sep="")[1:5000]
|
|
| 49 | 32 |
|
| 50 | 33 |
## vars to process |
| 51 | 34 |
vars=as.data.frame(matrix(c( |
| 52 | 35 |
"Cloud_Mask", "CM", "NN", 1, |
| 53 |
# "Sensor_Azimuth", "ZA", "CUBIC", 1, |
|
| 54 | 36 |
"Sensor_Zenith", "SZ", "CUBIC", 1), |
| 55 | 37 |
byrow=T,ncol=4,dimnames=list(1:2,c("variable","varid","method","band"))),stringsAsFactors=F)
|
| 56 | 38 |
|
| ... | ... | |
| 59 | 41 |
gpp = SpatialPolygons(list(Polygons(list(Polygon(gbb)),1))) |
| 60 | 42 |
proj4string(gpp)=projection(glb) |
| 61 | 43 |
|
| 62 |
outdir="~/acrobates/adamw/projects/interp/data/modis/mod35/processpath/gridded/"
|
|
| 44 |
outdir="/gridded/" |
|
| 63 | 45 |
|
| 64 | 46 |
swtif<-function(file,var){
|
| 65 | 47 |
outfile=paste(tempdir(),"/",var$varid,"_",basename(file),sep="") #gridded path |
| ... | ... | |
| 78 | 60 |
RESAMPLING_TYPE =",var$method," |
| 79 | 61 |
OUTPUT_PROJECTION_TYPE = GEO |
| 80 | 62 |
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 | 63 |
ELLIPSOID_CODE = WGS84 |
| 84 | 64 |
OUTPUT_TYPE = HDFEOS |
| 85 | 65 |
OUTPUT_FILENAME= ",outfile," |
| ... | ... | |
| 91 | 71 |
## now run the swath2grid tool |
| 92 | 72 |
## write the gridded file |
| 93 | 73 |
print(paste("Starting",file))
|
| 94 |
system(paste("sudo ",stitch," -p ",tempdir(),"/",basename(file),"_MODparms.txt -d -log /dev/null ",sep=""),intern=F)#,ignore.stderr=F)
|
|
| 74 |
system(paste("sudo ",stitch," -p ",tempdir(),"/",basename(file),"_MODparms.txt -d -log /dev/null -tmpLatLondir ",tempdir(),sep=""),intern=F)#,ignore.stderr=F)
|
|
| 95 | 75 |
print(paste("Finished processing variable",var$variable,"from ",basename(file),"to",outfile))
|
| 96 | 76 |
} |
| 97 | 77 |
|
| ... | ... | |
| 138 | 118 |
else return(0) |
| 139 | 119 |
})) |
| 140 | 120 |
|
| 121 |
## report on any bad files |
|
| 141 | 122 |
table(check) |
| 142 | 123 |
|
| 124 |
## remove any fail the check |
|
| 143 | 125 |
file.remove(gfiles[check==0]) |
| 126 |
gfiles=gfiles[check==1] |
|
| 144 | 127 |
|
| 145 |
## use new gdal |
|
| 146 |
#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 |
|
|
| 150 |
|
|
| 151 |
### Merge them into a geotiff |
|
| 152 |
#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 |
|
|
| 155 |
|
|
| 156 |
## try with pktools |
|
| 157 |
## global |
|
| 158 |
#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 |
#bb="-ulx -180 -uly 90 -lrx 180 -lry -90" |
|
| 160 |
#bb="-ulx -180 -uly 90 -lrx 170 -lry 80" |
|
| 161 |
#bb="-ulx -72 -uly 11 -lrx -59 -lry -1" |
|
| 162 |
|
|
| 163 |
|
|
| 164 |
#expand.grid(x=seq(-180,170,by=10),y=seq(-90,80)) |
|
| 165 |
#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 |
|
|
| 168 |
# bounding box? |
|
| 169 | 128 |
|
| 170 | 129 |
########### |
| 171 |
### Use GRASS to import all the tifs and calculat the mode |
|
| 130 |
### Use GRASS to import all the tifs and calculate the mode
|
|
| 172 | 131 |
## make temporary working directory |
| 173 | 132 |
tf=paste(tempdir(),"/grass", Sys.getpid(),"/", sep="") #temporar |
| 174 | 133 |
if(!file.exists(tf)) dir.create(tf) |
| ... | ... | |
| 190 | 149 |
## read in all tifs |
| 191 | 150 |
for(f in gfiles[!imported]) {
|
| 192 | 151 |
print(f) |
| 193 |
execGRASS("r.in.gdal",input=f,output=basename(f),flags="o")
|
|
| 152 |
execGRASS("r.external",input=f,output=basename(f),flags=c("overwrite"))
|
|
| 194 | 153 |
} |
| 195 | 154 |
|
| 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"))
|
|
| 155 |
## calculate mode in chunks. This first bins several individual swaths together into more or less complete global coverages taking the mode of each chunk |
|
| 156 |
nbreaks=100 |
|
| 157 |
bins=cut(1:5000,nbreaks) |
|
| 158 |
ts=system("g.mlist type=rast pattern=MOD*.tif",intern=T) #files to process
|
|
| 202 | 159 |
|
| 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"))
|
|
| 160 |
for(i in 1:nbreaks) #loop over breaks |
|
| 161 |
execGRASS("r.series",input=paste(ts[bins==levels(bins)[i]],sep="",collapse=","),output=paste("path",i,sep="_"),method="mode",range=c(0,5),flags=c("verbose","overwrite"),Sys_wait=T)
|
|
| 162 |
|
|
| 163 |
## Get mode of each chunk |
|
| 164 |
execGRASS("r.series",input=paste(system("g.mlist type=rast pattern=path*",intern=T),sep="",collapse=","),output="MOD35_path",method="mode",range=c(0,5),flags=c("verbose","overwrite"))
|
|
| 165 |
|
|
| 166 |
## fill in missing data (due to gridding artifacts) very near poles with water (north) and land (south) |
|
| 167 |
system("r.mapcalc \"MOD35_patha=if(isnull(MOD35_path)&y()>-84.31,0,MOD35_path)\"")
|
|
| 168 |
system("r.mapcalc \"MOD35_pathb=if(isnull(MOD35_patha)&y()<-84.31,3,MOD35_patha)\"")
|
|
| 205 | 169 |
|
| 206 | 170 |
## add colors |
| 207 |
execGRASS("r.colors",map="MOD35_path",rules="MOD35_path_grasscolors.txt")
|
|
| 171 |
execGRASS("r.colors",map="MOD35_pathb",rules="MOD35_path_grasscolors.txt")
|
|
| 172 |
|
|
| 208 | 173 |
## 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")
|
|
| 174 |
execGRASS("r.out.gdal",input="MOD35_pathb",output=paste(getwd(),"/C5MOD35_ProcessPath.tif",sep=""),type="Byte",createopt="COMPRESS=LZW,LEVEL=9,PREDICTOR=2")
|
|
| 175 |
|
|
| 176 |
## update metadata |
|
| 177 |
tags=c("TIFFTAG_IMAGEDESCRIPTION='Collection 5 MOD35 Processing Path (0=Water,1=Coast,2=Desert,3=Land)'",
|
|
| 178 |
"TIFFTAG_DOCUMENTNAME='Collection 5 MOD35 Processing Path'", |
|
| 179 |
"TIFFTAG_DATETIME='20130901'", |
|
| 180 |
"TIFFTAG_ARTIST='Adam M. Wilson (adam.wilson@yale.edu)'") |
|
| 181 |
system(paste("/usr/local/src/gdal-1.10.0/swig/python/scripts/gdal_edit.py ",getwd(),"/C5MOD35_ProcessPath.tif ",paste("-mo ",tags,sep="",collapse=" "),sep=""))
|
|
| 210 | 182 |
|
| 211 | 183 |
### delete the temporary files |
| 212 | 184 |
unlink_.gislock() |
| 213 | 185 |
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") |
|
| 229 |
|
|
Also available in: Unified diff
Updated C5 MOD35 Processing Path code and C5MOD35_Evaluation in preparation for publication