| 1 |
d39ab57e
|
Adam M. Wilson
|
## Short script to test swtif gridding of various versions
|
| 2 |
|
|
setwd(tempdir())
|
| 3 |
|
|
|
| 4 |
|
|
|
| 5 |
|
|
library(sp)
|
| 6 |
|
|
library(raster)
|
| 7 |
|
|
|
| 8 |
|
|
system("wget ftp://ladsweb.nascom.nasa.gov/allData/6/MOD35_L2/2009/029/MOD35_L2.A2009029.0500.006.2012245113542.hdf")
|
| 9 |
|
|
system("wget ftp://ladsweb.nascom.nasa.gov/allData/6/MOD35_L2/2009/029/MOD35_L2.A2009029.0320.006.2012245113606.hdf")
|
| 10 |
|
|
|
| 11 |
|
|
files=list.files(pattern="hdf")
|
| 12 |
|
|
|
| 13 |
|
|
swtifpath="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"
|
| 14 |
|
|
swtifpath="sudo MRTDATADIR=\"/usr/local/heg/2.12b/data\" PGSHOME=/usr/local/heg/2.12b/TOOLKIT_MTD PWD=/home/adamw /usr/local/heg/2.12b/bin/swtif"
|
| 15 |
|
|
|
| 16 |
|
|
## global bounding box
|
| 17 |
|
|
gbb=cbind(c(-180,-180,180,180,-180),c(-90,90,90,-90,-90))
|
| 18 |
|
|
gpp = SpatialPolygons(list(Polygons(list(Polygon(gbb)),1)))
|
| 19 |
|
|
proj4string(gpp)="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs "
|
| 20 |
|
|
|
| 21 |
|
|
vars=as.data.frame(matrix(c(
|
| 22 |
|
|
# "Cloud_Mask", "CM", "NN", 1,
|
| 23 |
|
|
# "Quality_Assurance", "QA", "NN", 1,
|
| 24 |
|
|
# "Solar_Zenith", "SolZen", "NN", 1,
|
| 25 |
|
|
"Sensor_Zenith", "SenZen", "CUBIC", 1
|
| 26 |
|
|
),
|
| 27 |
|
|
byrow=T,ncol=4,dimnames=list(1:1,c("variable","varid","method","band"))),stringsAsFactors=F)
|
| 28 |
|
|
|
| 29 |
|
|
|
| 30 |
|
|
## define function that grids swaths
|
| 31 |
|
|
swtif<-function(file,var){
|
| 32 |
|
|
outfile=paste(tempdir(),"/",var$varid,"_",basename(file),sep="") #gridded path
|
| 33 |
|
|
## First write the parameter file (careful, heg is very finicky!)
|
| 34 |
|
|
hdr=paste("NUM_RUNS = 1")
|
| 35 |
|
|
grp=paste("
|
| 36 |
|
|
BEGIN
|
| 37 |
|
|
INPUT_FILENAME=",file,"
|
| 38 |
|
|
OBJECT_NAME=mod35
|
| 39 |
|
|
FIELD_NAME=",var$variable,"|
|
| 40 |
|
|
BAND_NUMBER = ",var$band,"
|
| 41 |
|
|
OUTPUT_PIXEL_SIZE_X=926.6
|
| 42 |
|
|
OUTPUT_PIXEL_SIZE_Y=926.6
|
| 43 |
|
|
# MODIS 1km Resolution
|
| 44 |
|
|
SPATIAL_SUBSET_UL_CORNER = ( ",bbox(gpp)[2,2]," ",bbox(gpp)[1,1]," )
|
| 45 |
|
|
SPATIAL_SUBSET_LR_CORNER = ( ",bbox(gpp)[2,1]," ",bbox(gpp)[1,2]," )
|
| 46 |
|
|
RESAMPLING_TYPE =",var$method,"
|
| 47 |
|
|
OUTPUT_PROJECTION_TYPE = SIN
|
| 48 |
|
|
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 )
|
| 49 |
|
|
# projection parameters from http://landweb.nascom.nasa.gov/cgi-bin/QA_WWW/newPage.cgi?fileName=sn_gctp
|
| 50 |
|
|
ELLIPSOID_CODE = WGS84
|
| 51 |
|
|
OUTPUT_TYPE = HDFEOS
|
| 52 |
|
|
OUTPUT_FILENAME = ",outfile,"
|
| 53 |
|
|
END
|
| 54 |
|
|
",sep="")
|
| 55 |
|
|
## write it to a file
|
| 56 |
|
|
cat(c(hdr,grp) , file=paste(tempdir(),"/",basename(file),"_MODparms.txt",sep=""))
|
| 57 |
|
|
## now run the swath2grid tool
|
| 58 |
|
|
## write the gridded file
|
| 59 |
|
|
system(paste(swtifpath," -p ",tempdir(),"/",basename(file),"_MODparms.txt -d -tmpLatLondir ",tempdir(),sep=""),intern=F,ignore.stderr=F)
|
| 60 |
|
|
print(paste("Finished processing variable",var$variable,"from ",basename(file),"to",outfile))
|
| 61 |
|
|
}
|
| 62 |
|
|
|
| 63 |
|
|
swtifpath="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"
|
| 64 |
|
|
swtifpath="sudo MRTDATADIR=\"/usr/local/heg/2.12b/data\" PGSHOME=/usr/local/heg/2.12b/TOOLKIT_MTD PWD=/home/adamw /usr/local/heg/2.12b/bin/swtif"
|
| 65 |
|
|
|
| 66 |
|
|
## run it
|
| 67 |
|
|
lapply(1:nrow(vars),function(i) swtif(files[1],vars[i,]))
|
| 68 |
|
|
lapply(1:nrow(vars),function(i) swtif(files[2],vars[i,]))
|
| 69 |
|
|
|
| 70 |
|
|
### make a plot to compare versions
|
| 71 |
|
|
library(rasterVis)
|
| 72 |
|
|
|
| 73 |
|
|
d=stack(list.files(pattern="SenZen.*hdf$"))
|
| 74 |
|
|
|
| 75 |
|
|
png(file="swtifCompareVersions.png",width=2000,height=1000)
|
| 76 |
|
|
levelplot(d,at=seq(0,100,len=100),col.regions=rainbow(100))
|
| 77 |
|
|
dif=d[[1]]-d[[2]]
|
| 78 |
|
|
levelplot(dif,at=seq(min(dif,na.rm=T),100,len=100),col.regions=rainbow(100))
|
| 79 |
|
|
dev.off()
|
