[R-sig-Geo] Pixel reference location?

[Robert J. Hijmans]

Lyndon,

This is where I see that you use crop

to.crop <- raster(out.name)
cropped <- crop(to.crop, crop.img, filename = out.name, overwrite = TRUE)

You show crop.img and raster(out.name) _after_ this is run (so I do
not know what "to.crop" looks like; and that matters; I find what you
do confusing because your output is overwriting your input file
out.name). Anyway, you are surprised that the extent of the output and
crop.img not the same, right? I am guessing that this is as it should
be. Crop takes a subset (columns and rows) of an existing raster,
given an 'Extent' object you provide (in your case the Extent object
is extracted form the RasterLayer). If you provide an Exent that does
not match the coordinates of the rows and columns of a RasterLayer, it
is adjusted to the nearest rows/columns, as you can not have half
cells. So I am assuming that raster(out.name) (after the function) has
the same extent as "to.crop", and that it should be, because crop.img
is within half a cell.

Given all this, it also seems that what you really want is something
much simpler, rather than crop, just set the Extent to what it should
be:

extent(to.crop) <- extent(crop.img)

Other cases might use resample, but it seems that here all you need to
do is just set the values of the object's extent to what they should
be.

> was under the impression that raster refers to the upper left of the
> upper left pixel in defining the spatial extent, but I am getting

That is correct. raster uses the upper left, lower left, upper right,
and lower right of the entire structure (the extreme pixels if you
like).

Best, Robert

On Fri, Nov 12, 2010 at 9:12 AM, Lyndon Estes <lestes at princeton.edu> wrote:
> Dear R-Sig-Geo,
>
> I have a question regarding pixel referencing with the raster package,
> from which I am using the crop function to resize MODIS imagery.  I
> was under the impression that raster refers to the upper left of the
> upper left pixel in defining the spatial extent, but I am getting
> results that seem to suggest that the center of the pixel is being
> used.
>
> What I am doing: I have written a function that calls the MODIS
> HEGtools to read-in, stitch, and reproject the image into an Albers
> Equal Area Conic projection.  The HEGtool has the capacity to subset
> the image, but the resulting subset image is not well rectified with
> the result one gets when stitching two full size images together
> (probably because of how I have defined projection parameters), so I
> am instead adding an extra step that resizes the image with the crop
> function, using a subset of the image that I created using ENVI.
>
> Here are the values for the crop image:
>
> crop.img
> class       : RasterLayer
> filename    : MOD250m_extenttemplate
> nrow        : 6116
> ncol        : 4200
> ncell       : 25687200
> min value   : ?
> max value   : ?
> projection  : +proj=aea +lat_1=-18 +lat_2=-32 +lat_0=-30 +lon_0=24
> +x_0=0 +y_0=0 +ellps=clrk66 +datum=NAD27 +units=m +no_defs
> +nadgrids=@conus, at alaska, at ntv2_0.gsb, at ntv1_can.dat
> xmin        : -120863.9
> xmax        : 852092.8
> ymin        : -497022.6
> ymax        : 919787.6
> xres        : 231.6564
> yres        : 231.6564
>
> (Note the seemingly bizarre projection parameters, which is another
> story, and is created by the HEGtool, which is supposed to use WGS4,
> but instead seems to spit out NAD27).
>
> Now, after the crop function is run, my resulting cropped MODIS image
> ends up with the following parameters:
>
> raster(out.name)
> class       : RasterLayer
> filename    : MOD13Q1.A2005225.h20v11.mosaic4.NDVI.tif
> nrow        : 6116
> ncol        : 4200
> ncell       : 25687200
> min value   : ?
> max value   : ?
> projection  : +proj=aea +lat_1=-18 +lat_2=-32 +lat_0=-30 +lon_0=24
> +x_0=0 +y_0=0 +ellps=clrk66 +datum=NAD27 +units=m +no_defs
> +nadgrids=@conus, at alaska, at ntv2_0.gsb, at ntv1_can.dat
> xmin        : -120748.1
> xmax        : 852208.6
> ymin        : -496906.8
> ymax        : 919903.5
> xres        : 231.6564
> yres        : 231.6564
>
> Looking at the xmins of the two images, -120863.9 - -120748.1 =
> -115.8, which is equal to half a pixel.  Same for the ymaxs.
>
> So, is this an issue of where in the pixel raster assigns the UL
> reference coordinate, and, if so, is there a way to specify how the
> reference is being done?  If not, am I doing something else wrong?
>
> Thanks in advance for any advice (function code follows).
>
> Cheers, Lyndon
>
>
> ModStitch <- function(modpath, hegpath, modnames, cropimage) {
> # 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
>
>  # 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), ".mosaic4.",
> 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
>    write(paste("SPATIAL_SUBSET_UL_CORNER = ", "( ",
>                substr(mod.list[[1]][grep("GRID_UL_CORNER_LATLON=",
> mod.list[[1]])], 23, 48), " )", sep = ""),
>          filename, append = TRUE)
>
>    # Pulls coordinates for lower right from last header file
>    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)
>
>    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
>
>    # 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
>  }  # Close loop
> }  # End function
>
> ModStitch(modpath = path.ndvi.dat, hegpath = hegpath, modnames =
> modnames, cropimage = envi.img)
>
>
>
>
>
>
>
>
>
>
> --
> 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
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>