| 1 |
1 |
### this script generates empty tifs that align with the MODIS land tiles, but are available globally
|
| 2 |
|
|
| 3 |
|
|
| 4 |
|
## get MODIS gring coordingates
|
| 5 |
|
tb=read.table("http://modis-land.gsfc.nasa.gov/pdf/sn_gring_10deg.txt",skip=7,nrows=648,header=F,sep="")
|
| 6 |
|
colnames(tb)=c("iv","ih","ll_lon","ll_lat","ul_lon","ul_lat","ur_lon","ur_lat","lr_lon","lr_lat")
|
| 7 |
|
tb$tile=paste("h",sprintf("%02d",tb$ih),"v",sprintf("%02d",tb$iv),sep="")
|
| 8 |
|
## drop null tiles
|
| 9 |
|
tb=tb[!is.na(tb$ll_lon)&tb$ll_lon!=-999,]
|
| 10 |
|
## transform to polygons
|
| 11 |
|
getgrid=function(i){
|
| 12 |
|
print(i)
|
| 13 |
|
x=tb[i,]
|
| 14 |
|
gcols=matrix(as.numeric(c(x[3:10],x[3:4])),ncol=2,byrow=T)
|
| 15 |
|
gpol= Polygons(list(Polygon(gcols)),tb[i,"tile"])
|
| 16 |
|
return(gpol)
|
|
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){
|
|
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 |
print(paste("Making a raster grid for tile",tile))
|
|
29 |
## Earth Width (m)
|
|
30 |
ew=20015109.354*2
|
|
31 |
## Tile width or height = earth width / 36 = (20015109.354 + 20015109.354) / 36 = 1111950.5196666666 m
|
|
32 |
tw=ew/36
|
|
33 |
## number of cells (1200 for 1km data)
|
|
34 |
nc=1200
|
|
35 |
##1 km Cell size = tile width / cells
|
|
36 |
cs=tw/nc
|
|
37 |
##X coordinate of tile lower-left corner = -20015109.354 + horizontal tile number * tile width
|
|
38 |
llx= -(ew/2) + (h * tw)
|
|
39 |
lly= -(ew/4) + ((17-v) * tw)
|
|
40 |
## projection
|
|
41 |
proj="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs"
|
|
42 |
## make extent
|
|
43 |
ext=extent(llx,llx+(nc*cs),lly,lly+(nc*cs))
|
|
44 |
grid=raster(ext,nrows=nc,ncol=nc,crs=proj)
|
|
45 |
names(grid)=tile
|
|
46 |
return(grid)
|
| 17 |
47 |
}
|
| 18 |
48 |
|
| 19 |
|
modgrid=SpatialPolygons(lapply(1:nrow(tb),getgrid))
|
| 20 |
|
proj4string(modgrid)=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
|
| 21 |
49 |
|
| 22 |
|
## transform to sinusoidal
|
| 23 |
|
modgrids=spTransform(modgrid,CRS("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs "))
|
|
50 |
gs=expand.grid(h=0:35,v=0:17)
|
|
51 |
gs$tile=paste("h",sprintf("%02d",gs$h),"v",sprintf("%02d",gs$v),sep="")
|
|
52 |
|
|
53 |
modtiles=lapply(1:nrow(gs),function(i) maketile(gs$tile[i]))
|
|
54 |
names(modtiles)=gs$tile
|
| 24 |
55 |
|
| 25 |
56 |
## two tiles for testing
|
| 26 |
57 |
tiles=c("h22v02","h03v11")
|
| ... | ... | |
| 28 |
59 |
## path to MOD11A1 file for this tile to align grid/extent
|
| 29 |
60 |
gridfile=list.files("/nobackupp4/datapool/modis/MOD11A2.005/2009.01.01",pattern=paste(tiles[1],".*[.]hdf$",sep=""),recursive=T,full=T)[1]
|
| 30 |
61 |
td=raster(paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_8Day_1km_LST:LST_Day_1km",sep=""))
|
| 31 |
|
#projection(td)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs +datum=WGS84 +ellps=WGS84 "
|
| 32 |
|
|
| 33 |
|
gridfile2=list.files("/nobackupp4/datapool/modis/MOD11A2.005/2009.01.01",pattern=paste(tiles[2],".*[.]hdf$",sep=""),recursive=T,full=T)[1]
|
| 34 |
|
td2=raster(paste("HDF4_EOS:EOS_GRID:\"",gridfile2,"\":MODIS_Grid_8Day_1km_LST:LST_Day_1km",sep=""))
|
| 35 |
|
|
| 36 |
|
|
Finished function to make raster extent for any modis tile