[R] Hdf files Download

David Winsemius dwinsemius at comcast.net
Tue Oct 22 21:43:00 CEST 2013 

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:test@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 = ( 500000.0 -1800000.0 )', filename, append=TRUE)
  write('SPATIAL_SUBSET_LR_CORNER = ( 1600000.0 -2700000.0 )', filename, append=TRUE)
  write('  ', filename, append=TRUE)
  write(paste('OUTPUT_FILENAME = ', workd, 'tmp', dirname1, '.tif', sep=""), filename, append=TRUE)
  write('  ', filename, append=TRUE)
  write('RESAMPLING_TYPE = NEAREST_NEIGHBOR', filename, append=TRUE)
  write('  ', filename, append=TRUE)
  write('OUTPUT_PROJECTION_TYPE = PS', filename, append=TRUE)
  write('  ', filename, append=TRUE)
  write('OUTPUT_PROJECTION_PARAMETERS = ( ', filename, append=TRUE)
  write(' 6378137.0 6356752.314245179 0.0', filename, append=TRUE)
  write(' 0.0 0.0 71.0', filename, append=TRUE)
  write(' 0.0 0.0 0.0', filename, append=TRUE)
  write(' 0.0 0.0 0.0', filename, append=TRUE)
  write(' 0.0 0.0 0.0 )', filename, append=TRUE)
  write('  ', filename, append=TRUE)
  write('DATUM = WGS84', filename, append=TRUE)
  write('  ', filename, append=TRUE)
  write('OUTPUT_PIXEL_SIZE = 250', filename, append=TRUE)
  write('  ', filename, append=TRUE)
  close(filename)
  
  # Mosaic the images to get the whole area:
  shell(cmd=paste(MRT, 'resample -p ', MRT, 'mrt', dirname1, '.prm', sep=""))
  # delete all hdf files!
  unlink(paste(getwd(), '/', BLOCK1, sep=""))
  unlink(paste(getwd(), '/', BLOCK2, sep=""))
  #unlink(paste(getwd(), '/', BLOCK3, sep=""))
  #unlink(paste(getwd(), '/', BLOCK4, sep=""))
}
# end of script!

-- 
David.

> -Roy
> 
> On Oct 22, 2013, at 9:40 AM, "Tereza Smejkalova" <terka424 at gmail.com> wrote:
> 
>> I am triing one more time since my first message was rejected by the filter:
>> 
>> 
>> 
>> Dear all, 
>> I need to download and process large amounts of MODIS surface reflectance
>> imagery. I have adapted a script written by T. Hengl (
>> <http://www.spatial-analyst.net/wiki/?title=Download_and_resampling_of_MODIS
>> _images>
>> http://www.spatial-analyst.net/wiki/?title=Download_and_resampling_of_MODIS_
>> images) for download from the new HTTP (since FTP is no longer available).
>> The downloaded .HDF files however have no headers and it is impossible to
>> work with them. If I download directly from the webpage everything is fine.
>> Does anyone know what the problem might be, please? 
>> 
>> I am attaching my code
>> 
>> 
>> Please it is urgent. 

It only delays response if you don't "read the directions" for your toys that "may require assembly".

-- 
David Winsemius
Alameda, CA, USA

More information about the R-help mailing list