1
|
### this script generates empty tifs that align with the MODIS land tiles, but are available globally
|
2
|
library(raster)
|
3
|
|
4
|
######
|
5
|
### proceed by raster using information at http://code.env.duke.edu/projects/mget/wiki/SinusoidalMODIS
|
6
|
maketile=function(tile,type="raster",outdir="MODTILES"){
|
7
|
## check tile
|
8
|
## list of tile outside valid region
|
9
|
nodata=c('h00v00','h01v00','h02v00','h03v00','h04v00','h05v00','h06v00','h07v00','h08v00','h09v00','h10v00',
|
10
|
'h11v00','h12v00','h13v00','h22v00','h23v00','h24v00','h25v00','h26v00','h27v00','h28v00','h29v00','h30v00','h31v00',
|
11
|
'h32v00','h33v00','h34v00','h35v00','h00v01','h01v01','h02v01','h03v01','h04v01','h05v01','h06v01','h07v01','h08v01',
|
12
|
'h09v01','h10v01','h25v01','h26v01','h27v01','h28v01','h29v01','h30v01','h31v01','h32v01','h33v01','h34v01','h35v01',
|
13
|
'h00v02','h01v02','h02v02','h03v02','h04v02','h05v02','h06v02','h07v02','h08v02','h27v02','h28v02','h29v02',
|
14
|
'h30v02','h31v02','h32v02','h33v02','h34v02','h35v02','h00v03','h01v03','h02v03','h03v03','h04v03','h05v03','h30v03',
|
15
|
'h31v03','h32v03','h33v03','h34v03','h35v03','h00v04','h01v04','h02v04','h03v04','h32v04','h33v04','h34v04','h35v04',
|
16
|
'h00v05','h01v05','h34v05','h35v05','h00v06','h35v06','h00v11','h35v11','h00v12','h01v12','h34v12','h35v12','h00v13',
|
17
|
'h01v13','h02v13','h03v13','h32v13','h33v13','h34v13','h35v13','h00v14','h01v14','h02v14','h03v14','h04v14','h05v14',
|
18
|
'h30v14','h31v14','h32v14','h33v14','h34v14','h35v14','h00v15','h01v15','h02v15','h03v15','h04v15','h05v15','h06v15',
|
19
|
'h07v15','h08v15','h27v15','h28v15','h29v15','h30v15','h31v15','h32v15','h33v15','h34v15','h35v15','h00v16','h01v16',
|
20
|
'h02v16','h03v16','h04v16','h05v16','h06v16','h07v16','h08v16','h09v16','h10v16','h25v16','h26v16','h27v16','h28v16',
|
21
|
'h29v16','h30v16','h31v16','h32v16','h33v16','h34v16','h35v16','h00v17','h01v17','h02v17','h03v17','h04v17','h05v17',
|
22
|
'h06v17','h07v17','h08v17','h09v17','h10v17','h11v17','h12v17','h13v17','h22v17','h23v17','h24v17','h25v17','h26v17',
|
23
|
'h27v17','h28v17','h29v17','h30v17','h31v17','h32v17','h33v17','h34v17','h35v17')
|
24
|
if(tile%in%nodata) {print(paste(tile," not in region with data, skipping"));return(NULL)}
|
25
|
##
|
26
|
h=as.numeric(substr(tile,2,3))
|
27
|
v=as.numeric(substr(tile,5,6))
|
28
|
## Earth Width (m)
|
29
|
ew=20015109.354*2
|
30
|
## Tile width or height = earth width / 36 = (20015109.354 + 20015109.354) / 36 = 1111950.5196666666 m
|
31
|
tw=ew/36
|
32
|
## number of cells (1200 for 1km data)
|
33
|
nc=1200
|
34
|
##1 km Cell size = tile width / cells
|
35
|
cs=tw/nc
|
36
|
##X coordinate of tile lower-left corner = -20015109.354 + horizontal tile number * tile width
|
37
|
llx= -(ew/2) + (h * tw)
|
38
|
lly= -(ew/4) + ((17-v) * tw)
|
39
|
## projection
|
40
|
proj="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs"
|
41
|
## make extent
|
42
|
ext=extent(llx,llx+(nc*cs),lly,lly+(nc*cs))
|
43
|
grid=raster(ext,nrows=nc,ncol=nc,crs=proj)
|
44
|
names(grid)=tile
|
45
|
if(type=="raster") {
|
46
|
writeRaster(grid,paste(outdir,"/",tile,".tif",sep=""),options=c("COMPRESS=LZW"),datatype="INT1U",overwrite=T)
|
47
|
return(grid)
|
48
|
}
|
49
|
if(type=="extent") {
|
50
|
grid2=data.frame(tile=tile,t(as.vector(extent(grid))))
|
51
|
names(grid2)=c("tile","xmin", "xmax", "ymin", "ymax")
|
52
|
return(grid2)
|
53
|
}
|
54
|
}
|
55
|
|
56
|
## run it and save the extents as an Rdata object for easy retrieval
|
57
|
gs=expand.grid(h=0:35,v=0:17)
|
58
|
gs$tile=paste("h",sprintf("%02d",gs$h),"v",sprintf("%02d",gs$v),sep="")
|
59
|
modtiles=lapply(1:nrow(gs),function(i) maketile(gs$tile[i]))
|
60
|
names(modtiles)=gs$tile
|
61
|
|
62
|
## function to confirm that this method results in identical (e same extent, number of rows and columns, projection, resolution, and origin) rasters as MODIS LST data
|
63
|
checkgrid<-function(tile){
|
64
|
print(tile)
|
65
|
## path to MOD11A1 file for this tile to align grid/extent
|
66
|
gridfile=list.files("/nobackupp4/datapool/modis/MOD11A2.005/2009.01.01",pattern=paste(tile,".*[.]hdf$",sep=""),recursive=T,full=T)[1]
|
67
|
if(!file.exists(gridfile)) return(NULL)
|
68
|
td=raster(paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_8Day_1km_LST:LST_Day_1km",sep=""))
|
69
|
td2=maketile(tile)
|
70
|
return(compareRaster(td,td2,res=T,orig=T,rotation=T))
|
71
|
}
|
72
|
|
73
|
testing=F
|
74
|
if(testing){
|
75
|
## make the comparison for every tile
|
76
|
dat=do.call(c,lapply(gs$tile,checkgrid))
|
77
|
## summarize the results
|
78
|
table(dat)
|
79
|
}
|
80
|
|
81
|
### make table of extents
|
82
|
modtiles=do.call(rbind,lapply(1:nrow(gs),function(i) maketile(gs$tile[i],type="extent")))
|
83
|
|
84
|
## get list of tiles with some land
|
85
|
land=read.table("http://modis-land.gsfc.nasa.gov/pdf/sn_bound_10deg.txt",skip=6,nrows=648,header=T)
|
86
|
land$tile=paste("h",sprintf("%02d",land$ih),"v",sprintf("%02d",land$iv),sep="")
|
87
|
land=land[land$lon_min!=-999,]
|
88
|
modtiles$land=modtiles$tile%in%land$tile
|
89
|
|
90
|
|
91
|
## write the table as text
|
92
|
write.csv(modtiles,"MODIS_Tiles.csv",row.names=F,quote=F)
|