Hi Robert, Many thanks for helping further with this.
To send you files I produced, I setup a shared folder for you on our webspace, which will have generated an email notification with a link in it. This should work for downloading. This contains the image that served as "to.crop" (MOD13Q1.A2005225.h20v11.mosaic.NDVI.tif) and "crop.img" (MOD250m_extenttemplate). The original MODIS images can be obtained here. f1 <- "ftp://e4ftl01.cr.usgs.gov/MOLT/MOD13Q1.005/2005.08.13/MOD13Q1.A2005225.h20v11.005.2008213005356.hdf"; f2 <- "ftp://e4ftl01.cr.usgs.gov/MOLT/MOD13Q1.005/2005.08.13/MOD13Q1.A2005225.h20v12.005.2008213055406.hdf"; download.file(f1, destfile = paste(testpath, BLOCK1, sep = "")) download.file(f2, destfile = paste(testpath, BLOCK2, sep = "")) If you want to try and recreate what I have done, the HEGtool for MODIS is here: # Download for HEGtool, should you want to use it http://newsroom.gsfc.nasa.gov/sdptoolkit/HEG/HEGDownload.html And then the code I used to process the data together with the HEGtool follows. # Conversion code #ModStitch <- function(hegpath, modnames, cropimage) { ModStitch <- function(hegpath, modnames) { # Uses the HEGtool from MODIS LPDAAC to read in, stitch, and reproject the MODIS NDVI and VI quality bands # References: http://spatial-analyst.net/book/MODIS_download_HR.R # Args: # modpath: Path to directory holding MODIS data # hegpath: Path to the HEGtool installation # modnames: Vector containing MODIS image names to process # ###crop image: A raster defining the extent to which the images should be resized # Returns: # Separate stitched MODIS NDVI and VI Quality bands, reprojected to Albers Equal Area conic (for South # Africa), saved as a GEOtiff library(raster) # Set tool and environment variables hegtool <- paste(hegpath, "/heg/bin/hegtool", sep = "") grd.stitch.tool <- paste(hegpath, "/heg/bin/subset_stitch_grid", sep = "") pgs.loc <- paste(hegpath, "/heg/TOOLKIT_MTD", sep = "") mrt.data.loc <- paste(hegpath, "/heg/data", sep = "") Sys.setenv(PGSHOME = pgs.loc) # Environment variable for HEGtool Sys.setenv(MRTDATADIR = mrt.data.loc) # Environment variable for HEGtool # Set path to MODIS data directory #setwd(modpath) # And process the specified MODIS images using hegtool to capture header file information mod.list <- vector("list", length(modnames)) # List to capture header file output # Call crop image #crop.img <- raster(cropimage) # Create and capture HegHdr.hdr files for (i in 1:length(mod.list)) { system(paste(hegtool, "-h", modnames[i]), intern = FALSE, wait = TRUE) mod.list[[i]] <- scan("HegHdr.hdr", what = "character", sep = "\n",) names(mod.list)[i] <- modnames[i] } # Then stitch and reproject both the NDVI layers and the error flag layers # Set up a couple of variables which.var <- c("250m 16 days NDVI|", "250m 16 days VI Quality|") which.abb <- c("NDVI", "qual") for(k in 1:length(which.var)) { # Loop through NDVI, then VI # First write the stitch parameters file # Stitch parameter file name pr.name <- paste("stitch_", k, ".prm", sep = "") # 1 = stitch for NDVI, 2 for VI quality out.name <- paste(substr(modnames[1], 1, 23), ".mosaic.", which.abb[k],".tif", sep = "") # Open connection to write stitch file filename <- file(pr.name, open = "wt") write("", filename, append = TRUE) write("NUM_RUNS = 1", filename, append = TRUE) # Not sure if this needs to be a variable write("", filename, append = TRUE) write("BEGIN", filename, append = TRUE) write("NUMBER_INPUTFILES = 2",filename, append = TRUE) # Alter this if more than 1 file is to be stitched write(paste("INPUT_FILENAMES = ", modnames[1], "|", modnames[2], sep = ""), filename, append = TRUE) write("OBJECT_NAME = MODIS_Grid_16DAY_250m_500m_VI|", filename, append = TRUE) write(paste("FIELD_NAME = ", which.var[k], sep = ""), filename, append = TRUE) write("BAND_NUMBER = 1", filename, append = TRUE) # Pulls coordinates for upper left from first header file (turn on if want the whole image) write(paste("SPATIAL_SUBSET_UL_CORNER = ", "( ", substr(mod.list[[1]][grep("GRID_UL_CORNER_LATLON=", mod.list[[1]])], 23, 48), " )", sep = ""), filename, append = TRUE) # Specified upper left corner pair (UL -21.74814913, 22.82538048) #write("SPATIAL_SUBSET_UL_CORNER = ( -21.74814913 22.82538048 )", filename, append = TRUE) # Pulls coordinates for lower right from last header file (turn this on if you want entire image) write(paste("SPATIAL_SUBSET_LR_CORNER = ", "( ", substr(mod.list[[length(mod.list)]][grep("GRID_LR_CORNER_LATLON=", mod.list[[length(mod.list)]])], 23, 48), " )", sep = ""), filename, append = TRUE) # Specified lower right coordinate pair (LR -34.22328403, 33.19532794) #write("SPATIAL_SUBSET_LR_CORNER = ( -34.22328403 33.19532794 )", filename, append = TRUE) write("OUTPUT_OBJECT_NAME = MODIS_Grid_16DAY_250m_500m_VI|", filename, append = TRUE) write(paste("OUTGRID_X_PIXELSIZE = ", substr(mod.list[[1]][grep("GRID_PIXEL_SIZE=", mod.list[[1]])], 17, 26), sep = ""), filename, append = TRUE) write(paste("OUTGRID_Y_PIXELSIZE = ", substr(mod.list[[1]][grep("GRID_PIXEL_SIZE=", mod.list[[1]])], 17, 26), sep = ""), filename, append = TRUE) write(paste("OUTPUT_FILENAME = ", out.name, sep = ""), filename, append = TRUE) write("SAVE_STITCHED_FILE = NO", filename, append = TRUE) write(paste("OUTPUT_STITCHED_FILENAME = sasquatch", k, ".hdf", sep = ""), filename, append = TRUE) write("OUTPUT_TYPE = GEO", filename, append = TRUE) write("RESAMPLING_TYPE = NN", filename, append = TRUE) write("OUTPUT_PROJECTION_TYPE = ALBERS", filename, append = TRUE) write("OUTPUT_PROJECTION_PARAMETERS = ( 0.0 0.0 -18.0 -32.0 24.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 )", filename, append = TRUE) write("END", filename, append = TRUE) write("", filename, append = TRUE) close(filename) # Then call the stitch tool and process it stitch.com <- paste(grd.stitch.tool, "-P", pr.name) system(stitch.com) # Works # Turning this off for now because it appears to be shifting the output image by referencing the center of # the pixel. # Then use Raster function to resize because subsetting with MODIS tool results in offset #to.crop <- raster(out.name) #cropped <- crop(to.crop, crop.img, filename = out.name, overwrite = TRUE) # Overwrites existing images # Try trim the image instead (doesn't work') #to.trim <- raster(out.name) #NAvalue(to.trim) <- -1 # Define NA values as less than 0 #print(paste("Starting raster trim for", which.var[k])) #out.name2 <- paste("trim", out.name, sep = "") #trim(to.trim, filename = out.name2, overwrite = TRUE, format = "GTiff") #print(paste("Finished raster trim for", which.var[k])) } # Close loop } # End function setwd(testpath) # Set to directory where MODIS imagery is. ModStitch(hegpath = hegpath, modnames = modnames) # Crop component turned off Many thanks, and I look forward to hearing what you find. Cheers, Lyndon -- Lyndon Estes Research Associate Woodrow Wilson School Princeton University +1-609-258-2392 (o) +1-609-258-6082 (f) +1-202-431-0496 (m) [email protected] _______________________________________________ R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo
