|
1
|
## 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()
|
