Revision 5f212fe8
Added by Adam Wilson over 11 years ago
| climate/procedures/MOD35_ExtractProcessPath.r | ||
|---|---|---|
| 39 | 39 |
|
| 40 | 40 |
#### Grid and mosaic the swath data |
| 41 | 41 |
|
| 42 |
stitch="sudo MRTDATADIR=\"/usr/local/heg/data\" PGSHOME=/usr/local/heg/TOOLKIT_MTD PWD=/home/adamw /usr/local/heg/bin/swtif" |
|
| 42 |
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" |
|
| 43 |
stitch="/usr/local/heg/2.12/bin/swtif" |
|
| 43 | 44 |
|
| 45 |
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" |
|
| 44 | 46 |
files=paste(getwd(),"/",list.files("swath",pattern="hdf$",full=T),sep="")
|
| 45 | 47 |
|
| 46 | 48 |
## vars to process |
| 47 | 49 |
vars=as.data.frame(matrix(c( |
| 48 |
"Cloud_Mask", "CM",
|
|
| 49 |
"Sensor_Azimuth", "ZA",
|
|
| 50 |
"Sensor_Zenith", "SZ"), |
|
| 51 |
byrow=T,ncol=2,dimnames=list(1:3,c("variable","varid"))),stringsAsFactors=F)
|
|
| 50 |
"Cloud_Mask", "CM", "NN", 1,
|
|
| 51 |
# "Sensor_Azimuth", "ZA", "CUBIC", 1,
|
|
| 52 |
"Sensor_Zenith", "SZ", "CUBIC", 1),
|
|
| 53 |
byrow=T,ncol=4,dimnames=list(1:2,c("variable","varid","method","band"))),stringsAsFactors=F)
|
|
| 52 | 54 |
|
| 53 | 55 |
## global bounding box |
| 54 |
gbb=cbind(c(-180,-180,180,180,-180),c(-85,85,85,-85,-85))
|
|
| 56 |
gbb=cbind(c(-180,-180,180,180,-180),c(-90,90,90,-90,-90))
|
|
| 55 | 57 |
gpp = SpatialPolygons(list(Polygons(list(Polygon(gbb)),1))) |
| 56 | 58 |
proj4string(gpp)=projection(glb) |
| 57 | 59 |
|
| 60 |
outdir="~/acrobates/adamw/projects/interp/data/modis/mod35/gridded/" |
|
| 58 | 61 |
|
| 59 |
getpath<- function(file){
|
|
| 60 |
setwd(tempdir()) |
|
| 61 |
bfile=sub(".hdf","",basename(file))
|
|
| 62 |
tempfile_path=paste(tempdir(),"/path_",basename(file),sep="") #gridded path |
|
| 63 |
tempfile_sz=paste(tempdir(),"/sz_",basename(file),sep="") # gridded sensor zenith |
|
| 64 |
tempfile2_path=paste(tempdir(),"/",bfile,".tif",sep="") #gridded/masked/processed path |
|
| 65 |
outfile=paste("~/acrobates/adamw/projects/interp/data/modis/mod35/gridded/",bfile,".tif",sep="") #final file
|
|
| 66 |
if(file.exists(outfile)) return(c(file,0)) |
|
| 67 |
## get bounding coordinates |
|
| 68 |
# glat=as.numeric(do.call(c,strsplit(sub("GRINGPOINTLATITUDE=","",system(paste("gdalinfo ",file," | grep GRINGPOINTLATITUDE"),intern=T)),split=",")))
|
|
| 69 |
# glon=as.numeric(do.call(c,strsplit(sub("GRINGPOINTLONGITUDE=","",system(paste("gdalinfo ",file," | grep GRINGPOINTLONGITUDE"),intern=T)),split=",")))
|
|
| 70 |
# bb=cbind(c(glon,glon[1]),c(glat,glat[1])) |
|
| 71 |
# pp = SpatialPolygons(list(Polygons(list(Polygon(bb)),1))) |
|
| 72 |
# proj4string(pp)=projection(glb) |
|
| 73 |
# ppc=gIntersection(pp,gpp) |
|
| 74 |
# ppc=gBuffer(ppc,width=0.3) #buffer a little to remove gaps between images |
|
| 75 |
ppc=gpp |
|
| 62 |
swtif<-function(file,var){
|
|
| 63 |
outfile=paste(tempdir(),"/",var$varid,"_",basename(file),sep="") #gridded path |
|
| 76 | 64 |
## First write the parameter file (careful, heg is very finicky!) |
| 77 |
hdr=paste("NUM_RUNS = ",length(vars$varid),"|MULTI_BAND_HDFEOS:",length(vars$varid),sep="")
|
|
| 65 |
hdr=paste("NUM_RUNS = 1")
|
|
| 78 | 66 |
grp=paste("
|
| 79 | 67 |
BEGIN |
| 80 | 68 |
INPUT_FILENAME=",file," |
| 81 | 69 |
OBJECT_NAME=mod35 |
| 82 |
FIELD_NAME=",vars$variable,"|
|
|
| 83 |
BAND_NUMBER = ",1:length(vars$varid),"
|
|
| 70 |
FIELD_NAME=",var$variable,"| |
|
| 71 |
BAND_NUMBER = ",var$band,"
|
|
| 84 | 72 |
OUTPUT_PIXEL_SIZE_X=0.008333333 |
| 85 | 73 |
OUTPUT_PIXEL_SIZE_Y=0.008333333 |
| 86 |
SPATIAL_SUBSET_UL_CORNER = ( ",bbox(ppc)[2,2]," ",bbox(ppc)[1,1]," ) |
|
| 87 |
SPATIAL_SUBSET_LR_CORNER = ( ",bbox(ppc)[2,1]," ",bbox(ppc)[1,2]," ) |
|
| 88 |
OUTPUT_OBJECT_NAME = mod35| |
|
| 89 |
RESAMPLING_TYPE =NN |
|
| 74 |
SPATIAL_SUBSET_UL_CORNER = ( ",bbox(gpp)[2,2]," ",bbox(gpp)[1,1]," ) |
|
| 75 |
SPATIAL_SUBSET_LR_CORNER = ( ",bbox(gpp)[2,1]," ",bbox(gpp)[1,2]," ) |
|
| 76 |
RESAMPLING_TYPE =",var$method," |
|
| 90 | 77 |
OUTPUT_PROJECTION_TYPE = GEO |
| 91 | 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 ) |
| 92 | 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 ) |
| 93 | 80 |
# projection parameters from http://landweb.nascom.nasa.gov/cgi-bin/QA_WWW/newPage.cgi?fileName=sn_gctp |
| 94 | 81 |
ELLIPSOID_CODE = WGS84 |
| 95 | 82 |
OUTPUT_TYPE = HDFEOS |
| 96 |
OUTPUT_FILENAME= ",tempfile_path,"
|
|
| 83 |
OUTPUT_FILENAME= ",outfile,"
|
|
| 97 | 84 |
END |
| 98 | 85 |
|
| 99 | 86 |
",sep="") |
| 100 |
|
|
| 101 |
## if any remnants from previous runs remain, delete them |
|
| 102 |
if(length(list.files(tempdir(),pattern=bfile)>0)) |
|
| 103 |
file.remove(list.files(tempdir(),pattern=bfile,full=T)) |
|
| 104 | 87 |
## write it to a file |
| 105 | 88 |
cat( c(hdr,grp) , file=paste(tempdir(),"/",basename(file),"_MODparms.txt",sep="")) |
| 106 | 89 |
## now run the swath2grid tool |
| 107 | 90 |
## write the gridded file |
| 108 | 91 |
print(paste("Starting",file))
|
| 109 |
system(paste("(",stitch," -p ",tempdir(),"/",basename(file),"_MODparms.txt -d)",sep=""),intern=F)#,ignore.stderr=F)
|
|
| 110 |
############## Now run the 5km summary |
|
| 111 |
## First write the parameter file (careful, heg is very finicky!) |
|
| 112 |
hdr=paste("NUM_RUNS = ",nrow(vars[-1,]),"|MULTI_BAND_HDFEOS:",nrow(vars[-1,]),sep="")
|
|
| 113 |
grp=paste("
|
|
| 114 |
BEGIN |
|
| 115 |
INPUT_FILENAME=",file," |
|
| 116 |
OBJECT_NAME=mod35 |
|
| 117 |
FIELD_NAME=",vars$variable[-1],"| |
|
| 118 |
BAND_NUMBER = ",1," |
|
| 119 |
OUTPUT_PIXEL_SIZE_X=0.008333333 |
|
| 120 |
#0.0416666 |
|
| 121 |
OUTPUT_PIXEL_SIZE_Y=0.008333333 |
|
| 122 |
#0.0416666 |
|
| 123 |
SPATIAL_SUBSET_UL_CORNER = ( ",bbox(ppc)[2,2]," ",bbox(ppc)[1,1]," ) |
|
| 124 |
SPATIAL_SUBSET_LR_CORNER = ( ",bbox(ppc)[2,1]," ",bbox(ppc)[1,2]," ) |
|
| 125 |
#OUTPUT_OBJECT_NAME = mod35| |
|
| 126 |
RESAMPLING_TYPE =NN |
|
| 127 |
OUTPUT_PROJECTION_TYPE = GEO |
|
| 128 |
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 ) |
|
| 129 |
#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 ) |
|
| 130 |
# projection parameters from http://landweb.nascom.nasa.gov/cgi-bin/QA_WWW/newPage.cgi?fileName=sn_gctp |
|
| 131 |
ELLIPSOID_CODE = WGS84 |
|
| 132 |
OUTPUT_TYPE = HDFEOS |
|
| 133 |
OUTPUT_FILENAME= ",tempfile_sz," |
|
| 134 |
END |
|
| 92 |
system(paste("",stitch," -p ",tempdir(),"/",basename(file),"_MODparms.txt -d -log /dev/null ",sep=""),intern=F)#,ignore.stderr=F)
|
|
| 93 |
print(paste("Finished processing variable",var$variable,"from ",basename(file),"to",outfile))
|
|
| 94 |
} |
|
| 135 | 95 |
|
| 136 |
",sep="") |
|
| 137 |
|
|
| 138 |
## write it to a file |
|
| 139 |
cat( c(hdr,grp) , file=paste(tempdir(),"/",basename(file),"_MODparms_angle.txt",sep="")) |
|
| 140 |
|
|
| 141 |
## now run the swath2grid tool |
|
| 142 |
## write the gridded file |
|
| 143 |
print(paste("Starting",file))
|
|
| 144 |
system(paste("(",stitch," -p ",tempdir(),"/",basename(file),"_MODparms_angle.txt -d)",sep=""),intern=F,ignore.stderr=F)
|
|
| 96 |
getpath<- function(file){
|
|
| 97 |
setwd(tempdir()) |
|
| 98 |
bfile=sub(".hdf","",basename(file))
|
|
| 99 |
tempfile2_path=paste(tempdir(),"/",bfile,".tif",sep="") #gridded/masked/processed path |
|
| 100 |
outfile=paste(outdir,"/",bfile,".tif",sep="") #final file |
|
| 101 |
if(file.exists(outfile)) return(c(file,0)) |
|
| 102 |
ppc=gpp |
|
| 103 |
####### |
|
| 104 |
## run swtif for each band |
|
| 105 |
lapply(1:nrow(vars),function(i) swtif(file,vars[i,])) |
|
| 145 | 106 |
####### import to R for processing |
| 146 |
if(!file.exists(tempfile_path)) {
|
|
| 107 |
|
|
| 108 |
if(!file.exists(paste(tempdir(),"/CM_",basename(file),sep=""))) {
|
|
| 147 | 109 |
file.remove(list.files(tempdir(),pattern=bfile,full=T)) |
| 148 | 110 |
return(c(file,0)) |
| 149 | 111 |
} |
| 150 | 112 |
## convert to land path |
| 151 |
d=raster(paste("HDF4_EOS:EOS_GRID:\"",tempfile_path,"\":mod35:Cloud_Mask_0",sep=""))
|
|
| 152 |
sz=raster(paste("HDF4_EOS:EOS_GRID:\"",tempfile_sz,"\":mod35:Sensor_Zenith",sep=""))
|
|
| 153 |
NAvalue(sz)=0 |
|
| 154 |
## resample sensor angles to 1km grid and mask paths with angles >=30 |
|
| 155 |
# sz2=resample(sz,d,method="ngb",file=sub("hdf","tif",tempfile_sz),overwrite=T)
|
|
| 156 |
getlc=function(x,y) {ifelse(y==0|y>40,NA,((x%/%2^6) %% 2^2))}
|
|
| 113 |
d=raster(paste(tempdir(),"/CM_",basename(file),sep="")) |
|
| 114 |
sz=raster(paste(tempdir(),"/SZ_",basename(file),sep="")) |
|
| 115 |
NAvalue(sz)=-9999 |
|
| 116 |
getlc=function(x,y) {ifelse(y==0|y>6000,NA,((x%/%2^6) %% 2^2))}
|
|
| 157 | 117 |
path= overlay(d,sz,fun=getlc,filename=tempfile2_path,options=c("COMPRESS=LZW", "LEVEL=9","PREDICTOR=2"),datatype="INT1U",overwrite=T)
|
| 158 | 118 |
### warp them to align all pixels |
| 159 | 119 |
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=""))
|
| ... | ... | |
| 163 | 123 |
} |
| 164 | 124 |
|
| 165 | 125 |
|
| 166 |
## establish sudo priveleges to run swtif
|
|
| 167 |
system("sudo ls"); mclapply(files,getpath,mc.cores=10)
|
|
| 126 |
### run it
|
|
| 127 |
mclapply(files[1:200],getpath,mc.cores=10)
|
|
| 168 | 128 |
|
| 169 | 129 |
## check gdal can read all of them |
| 170 |
gfiles=list.files("gridded",pattern="tif$",full=T)
|
|
| 130 |
gfiles=list.files(outdir,pattern="tif$",full=T)
|
|
| 171 | 131 |
length(gfiles) |
| 172 | 132 |
|
| 173 | 133 |
check=do.call(rbind,mclapply(gfiles,function(file){
|
| ... | ... | |
| 181 | 141 |
file.remove(gfiles[check==0]) |
| 182 | 142 |
|
| 183 | 143 |
## use new gdal |
| 184 |
system(paste("/usr/local/gdal-1.10.0/bin/gdalwarp -wm 900 -overwrite -co COMPRESS=LZW -co PREDICTOR=2 -multi -r mode gridded/*.tif MOD35_path_gdalwarp.tif"))
|
|
| 144 |
system(paste("/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.tif",sep=""))
|
|
| 185 | 145 |
|
| 186 | 146 |
|
| 187 | 147 |
### Merge them into a geotiff |
Also available in: Unified diff
Added script to test functionality of HEG tool which reveals that v2.12 segfaults when multiple bands are selected for processing. Also updated MOD35_ExtractProcessPath to process one-band-at-a-time instead of in batches.