[R-sig-Geo] Download MODIS images

[Tomislav Hengl]

> -----Original Message-----
> From: r-sig-geo-bounces at stat.math.ethz.ch [mailto:r-sig-geo-
> bounces at stat.math.ethz.ch] On Behalf Of Rosca Bogdan
> Sent: Saturday, July 03, 2010 7:04 PM
> To: r-sig-geo at stat.math.ethz.ch
> Subject: [R-sig-Geo] Download MODIS images
> 
> Hi!
> I am trying to work with T. Hengl's script downloaded from:
> http://spatial-analyst.net/book/Rscripts in order to get modis images for my
> country. I have only modify the script in order to sample data according to
> my country border (Romania, which is not so far away from Croatia) and
> everything run smooth until these lines:
> 
> # download the blocks in a loop (year 2006):
> for (i in 268:313){
> getlist <- strsplit(getURL(paste(MOD11A2, dates$dirname[[i]], "/", sep=""),
> .opts=curlOptions(ftplistonly=
> TRUE)), "\r\n")[[1]]
> BLOCK1 <- getlist[grep(getlist,
> pattern="MOD11A2.********.h18v04.*************.hdf")[1]]
> BLOCK2 <- getlist[grep(getlist,
> pattern="MOD11A2.********.h19v04.*************.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(MOD11A2a, dates$dirname[[i]], "/", BLOCK1,sep=""),
> destfile=paste(getwd(), "/", BLOCK1, sep=""), mode='wb', method='wget',
> quiet=T, cacheOK=FALSE)
> download.file(paste(MOD11A2, dates$dirname[[i]], "/", BLOCK2,sep=""),
> destfile=paste(getwd(), "/", BLOCK2,sep=""), mode='wb', method='wget',
> quiet=T, 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)
> close(mosaicname)
> # generate temporary mosaic:
> shell(cmd=paste(MRT, 'mrtmosaic -i ', MRT, 'TmpMosaic.prm -s "0 0 0 0 1 0 0
> 0 0 0 0 0" -o ', workd, 'TmpMosaic.hdf', sep=""))
> close(mosaicname)
> 
> # resample to UTM 33N:
> filename <- file(paste(MRT, "mrt", dirname1, ".prm", sep=""), open="wt")
> write(paste('INPUT_FILENAME = ', workd, 'TmpMosaic.hdf', sep=""), filename)
> 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(paste('SPATIAL_SUBSET_UL_CORNER = (', XUL, YUL, ')'), filename,
> append=TRUE)
> write(paste('SPATIAL_SUBSET_LR_CORNER = (', XLR, YLR, ')'), filename,
> append=TRUE)
> write('  ', filename, append=TRUE)
> write(paste('OUTPUT_FILENAME = ', workd, 'LST', 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 = UTM', filename, append=TRUE)
> write('  ', filename, append=TRUE)
> write('OUTPUT_PROJECTION_PARAMETERS = ( ', filename, append=TRUE)
> write(paste(lon.c, lat.c, '0.0'), filename, append=TRUE)
> write(' 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(' 0.0 0.0 0.0 )', filename, append=TRUE)
> write('  ', filename, append=TRUE)
> write('UTM_ZONE = 35', filename, append=TRUE)
> write('DATUM = WGS84', filename, append=TRUE)
> write('  ', filename, append=TRUE)
> write('OUTPUT_PIXEL_SIZE = 1000', 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=""))
> }
> 
> I keep getting an error like this:
> Error in grep(getlist, pattern =
> "MOD11A2.********.h18v04.*************.hdf") :
>  invalid regular expression 'MOD11A2.********.h18v04.*************.hdf'
> In addition: Warning message:
> In grep(getlist, pattern = "MOD11A2.********.h18v04.*************.hdf") :
>  regcomp error:  'Invalid use of repetition operators'

Modify to:

> BLOCK1 <- getlist[grep(getlist, pattern="MOD11A2.*.h18v04.*.hdf")[1]]

One more thing. Connection with the NASA's FTP often breaks (security setting of their server?) so that the derived object might finish being empty or incomplete.

T. Hengl


> Please help.
> Thank's,
> Rosca Bogdan
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo