[R-sig-Geo] Pixel reference location?

[Lyndon Estes]

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)
lestes at princeton.edu