On Oct 22, 2013, at 11:34 AM, Roy Mendelssohn - NOAA Federal wrote:
Your script didn't make it. There are restrictions on the mail list for
attached files. You might try putting the script directly into the email.
Thinking that the rejection might have be due to a Nabble posting I located it
there and here include it. I'm hoping to use the method as well.
#
# title : MODIS_download_HTTP.R
# purpose : Download of MODIS EVI images for the British Colombia;
# reference :
http://spatial-analyst.net/wiki/index.php?title=Download_and_resampling_of_MODIS_images
# producer : Prepared by T. Hengl (addapted by T. Smejkalova)
# last update : Southampton, UK, 18 October 2013.
# inputs: Coordinates of the area of interest; proj4 parameters; ftp
addresses etc.;
# outputs : a series of EVI images (geoTIFF) projected in the local
coordinate system;
# remarks 1 : To run this script, you need to obtain and install the MODIS
resampling tool from
[https://lpdaac.usgs.gov/lpdaac/tools/modis_reprojection_tool];
# remarks 2 : You should also obtain the WGET from
[http://users.ugent.be/~bpuype/wget/] --- simply copy the wget exe to windows
system folder;
# remarks 3 : make sure you disable your antivirus tools such as Norton or
McAfee otherwise it might block wget from running!
library(rgdal)
library(RCurl)
setInternet2(use = TRUE)
# Obtain the MODIS tool from:
http://lpdaac.usgs.gov/landdaac/tools/modis/index.asp
# setwd(E:/PineBeetleBC/MODIS)
# location of the MODIS 1 km monthly blocks:
MOD09GQ - http://e4ftl01.cr.usgs.gov/MOLT/MOD09GQ.005/;
MOD09GQa - http://anonymous:t...@e4ftl01.cr.usgs.gov/MOLT/MOD09GQ.005/;
product = MOD09GQ
year = 2000
# location of the mosaicing tool:
MRT - 'D:\\MRT\\MRT_Win\\bin\\'
workd - 'D:\\R_Modis\\'
setwd('D:\\R_Modis\\')
options(download.file.method=auto)
# get the list of directories (thanks to Barry Rowlingson):
items1 - strsplit(getURL(MOD09GQ), \n)[[1]]
# get the directory names and create a new data frame:
dates=data.frame(dirname=substr(items1[20:length(items1)],52,62))
# get the list of *.hdf files:
dates$BLOCK1 - rep(NA, length(dates))
dates$BLOCK2 - rep(NA, length(dates))
#dates$BLOCK3 - rep(NA, length(dates))
#dates$BLOCK4 - rep(NA, length(dates))
for (i in 1:1)#length(dates))# for each date per year
{
getlist - strsplit(getURL(paste(MOD09GQ, dates$dirname[i], sep=)),
)[[1]]
getlist=getlist[grep(product,getlist)]
getlist=getlist[grep(.hdf,getlist)]
filenames=substr(getlist,1,nchar(getlist[1])-3)
BLOCK1 - filenames[grep(filenames,pattern=MOD09GQ.*.h18v02.*.hdf)[1]]
BLOCK2 - filenames[grep(filenames,pattern=MOD09GQ.*.h19v02.*.hdf)[1]]
# write up the file names back to the dates.txt:
for(j in 2:3){
dates[i,j] - get(paste(BLOCK, j-1, sep=))
}
# Download all blocks from the list to a local drive:
#
while(!is.na(dates[i,2])!is.na(dates[i,3])!is.na(dates[i,4])!is.na(dates[i,5])!is.na(dates[i,6])!is.na(dates[i,7])!is.na(dates[i,8])!is.na(dates[i,9])!is.na(dates[i,10])){
download.file(paste(MOD09GQa, dates$dirname[i],
dates$BLOCK1[i],sep=),destfile=paste(getwd(), /, BLOCK1, sep=))
download.file(paste(MOD09GQ, dates$dirname[i],
dates$BLOCK2[i],sep=),destfile=paste(getwd(), /, BLOCK2, sep=),
cacheOK=FALSE, )
# remove . from the file name:
dirname1 - sub(sub(pattern=\\., replacement=_, dates$dirname[[i]]),
pattern=\\., replacement=_, dates$dirname[[i]])
# mosaic the blocks:
mosaicname = file(paste(MRT, TmpMosaic.prm, sep=), open=wt)
write(paste(workd, BLOCK1, sep=), mosaicname)
write(paste(workd, BLOCK2, sep=), mosaicname, append=T)
#write(paste(workd, BLOCK3, sep=), mosaicname, append=T)
#write(paste(workd, BLOCK4, sep=), mosaicname, append=T)
close(mosaicname)
# generate temporary mosaic:
shell(cmd=paste(MRT, 'mrtmosaic -i ', MRT, 'TmpMosaic.prm -s 0 1 0 0 0 0 0 0
0 0 0 -o ', workd, 'TmpMosaic.hdf', sep=))
# resample to epsg=3005:
filename = file(paste(MRT, mrt, dirname1, .prm, sep=), open=wt)
write(paste('INPUT_FILENAME = ', workd, 'TmpMosaic.hdf', sep=), filename)
# write(paste('INPUT_FILENAMES = ( ', workd, BLOCK1, ' ', workd, BLOCK2, ' ',
workd, BLOCK3, ' ', workd, BLOCK4, ' ', workd, BLOCK5, ' ', workd, BLOCK6, ' ',
workd, BLOCK7, ' ', workd, BLOCK8, ' ', workd, BLOCK9, ' )', sep=), filename)
# unfortunatelly does not work via command line :(
write(' ', filename, append=TRUE)
# write('SPECTRAL_SUBSET = ( 0 1 0 0 0 0 0 0 0 0 0 )', filename, append=TRUE)
write('SPECTRAL_SUBSET = ( 1 )', filename, append=TRUE)
write(' ', filename, append=TRUE)
write('SPATIAL_SUBSET_TYPE = OUTPUT_PROJ_COORDS', filename, append=TRUE)
write(' ', filename, append=TRUE)
write('SPATIAL_SUBSET_UL_CORNER = ( 50.0 -180.0 )', filename,
append=TRUE)
write('SPATIAL_SUBSET_LR_CORNER = ( 160.0 -270.0 )', filename,
append=TRUE)
write(' ', filename, append=TRUE)
write(paste('OUTPUT_FILENAME = ', workd,
