|
1
|
### summarize MOD35 data
|
|
2
|
setwd("/home/adamw/acrobates/adamw/projects/interp/")
|
|
3
|
|
|
4
|
library(sp)
|
|
5
|
library(spgrass6)
|
|
6
|
library(rgdal)
|
|
7
|
library(reshape)
|
|
8
|
library(ncdf4)
|
|
9
|
library(geosphere)
|
|
10
|
library(rgeos)
|
|
11
|
library(multicore)
|
|
12
|
library(raster)
|
|
13
|
library(lattice)
|
|
14
|
library(rgl)
|
|
15
|
library(hdf5)
|
|
16
|
library(rasterVis)
|
|
17
|
library(heR.Misc)
|
|
18
|
library(car)
|
|
19
|
|
|
20
|
X11.options(type="Xlib")
|
|
21
|
ncores=20 #number of threads to use
|
|
22
|
|
|
23
|
psin=CRS("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs")
|
|
24
|
|
|
25
|
tile="h11v08"
|
|
26
|
|
|
27
|
## get MODLAND tile information
|
|
28
|
tb=read.table("http://landweb.nascom.nasa.gov/developers/sn_tiles/sn_bound_10deg.txt",skip=6,nrows=648,header=T)
|
|
29
|
tb$tile=paste("h",sprintf("%02d",tb$ih),"v",sprintf("%02d",tb$iv),sep="")
|
|
30
|
save(tb,file="modlandTiles.Rdata")
|
|
31
|
|
|
32
|
tile="h11v08" #can move this to submit script if needed
|
|
33
|
#tile="h09v04" #oregon
|
|
34
|
|
|
35
|
tile_bb=tb[tb$tile==tile,] ## identify tile of interest
|
|
36
|
roi_ll=extent(tile_bb$lon_min,tile_bb$lon_max,tile_bb$lat_min,tile_bb$lat_max)
|
|
37
|
#roi=spTransform(roi,psin)
|
|
38
|
#roil=as(roi,"SpatialLines")
|
|
39
|
|
|
40
|
dmod06="data/modis/mod06/summary"
|
|
41
|
dmod35="data/modis/mod35/summary"
|
|
42
|
|
|
43
|
|
|
44
|
|
|
45
|
|
|
46
|
##################################################
|
|
47
|
### generate KML file for viewing in google earth
|
|
48
|
file=paste(dmod35,"/MOD35_",tile,"_ymonmean.nc",sep="")
|
|
49
|
m=1
|
|
50
|
|
|
51
|
cat("
|
|
52
|
0 green
|
|
53
|
50 grey
|
|
54
|
90 blue
|
|
55
|
100 red
|
|
56
|
nv 0 0 0 0
|
|
57
|
",file=paste(dmod35,"/mod35_colors.txt",sep=""))
|
|
58
|
|
|
59
|
system(paste("gdalinfo ",file,sep=""))
|
|
60
|
|
|
61
|
system(paste("gdalwarp -multi -r cubicspline -srcnodata 255 -dstnodata 255 -s_srs '+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs' -t_srs 'EPSG:4326' ",file, " MOD35_",tile,".tif",sep=""))
|
|
62
|
system(paste("gdaldem color-relief -b ",m," ",dmod35,"/MOD35_",tile,".tif mod35_colors.txt ",dmod35,"/MOD35_",tile,"_",m,".tif",sep=""))
|
|
63
|
system(paste("gdalinfo MOD35_",tile,"_",m,".tif",sep=""))
|
|
64
|
|
|
65
|
#system(paste("gdal_translate -b ",m," MOD35_",tile,".tif MOD35_",tile,"_",m,".tif",sep=""))
|
|
66
|
40075=256*(2^z)) #find zoom level
|
|
67
|
|
|
68
|
system(paste("gdal2tiles.py -z 1-8 -k ",dmod35,"/MOD35_",tile,"_",m,".tif",sep=""))
|
|
69
|
|
|
70
|
### Wind Direction
|
|
71
|
winddir="home/adamw/acrobates/adamw/projects/interp/data/ncep/"
|
|
72
|
dir.create(winddir)
|
|
73
|
system(paste("wget ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis.derived/surface/vwnd.mon.ltm.nc -P ",winddir))
|
