| 1 |
99b3508d
|
Adam M. Wilson
|
### 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))
|
