Pixel reference location?

Hi Robert,

Thanks for your response, and I am sorry I did not post the values for
"to.crop". Please see my responses below interspersed with your text.

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;

Here are to.crop's values:

to.crop
class ? ? ? : RasterLayer
filename ? ?MOD13Q1.A2005225.h20v11.mosaic4.NDVI.tif
nrow ? ? ? ?: 9872
ncol ? ? ? ?: 6976
ncell ? ? ? : 68867072
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 ? ? ? ?: -283139.2
xmax ? ? ? ?: 1332896
ymin ? ? ? ?: -1172417
ymax ? ? ? ?: 1114495
xres ? ? ? ?: 231.6564
yres ? ? ? ?: 231.6564

Interestingly, in both ENVI and ArcGIS the coordinate pair for the
upper left is listed as:

-283255.040, 1114610.641

In ArcGIS, I can tell that ArcGIS ascribes these coordinates to the
upper left corner of the pixel, where the o's below roughly describe
the outline of a pixel (in the hopes that the formatting makes it
through clearly):

xoooooo
ooooooo
ooooooo

where the x is located.

In ENVI, when I open the image in question, it tells me that the
Albers coordinate pair -283255.040, 1114610.641 is tied to the image
coordinates of 1.5, 1.5, which should be the center of the pixel.

ooooooo
oooxooo
ooooooo

I find what you

do confusing because your output is overwriting your input file
out.name).

Yes, I suppose that is bad practice. ?I was trying to cut down on disk
storage space, but I probably should write a new file and then use
file.remove to get rid of the larger image (which is about 1/3 ocean).

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.

I actually made the extent image out of "to.crop", by specifying
certain pixel coordinates in ENVI, but I now see that after making the
extent image, ENVI says that its image tie coordinates are 1, 1,
whereas to.crop's tie coordinates are 1.5, 1.5 (I had assumed that
they were at 1,1). ?I suspect that it is this problem that I am seeing
after cropping with R. However, please see my last question below:

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).

My question on this subject relates to the rough pixels I have
sketched above and the image tie coordinates (sorry if I was unclear).
?What I was curious to know is whether raster "thinks" that the
coordinate pair of

-283139.2, 1114495

is where the x is here (image coordinate 1,1, in case this doesn't
display well):

xoooooo
ooooooo
ooooooo

or here (image coordinate 1.5, 1.5)?

ooooooo
oooxooo
ooooooo

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)

Thanks, "extent" sounds much better, and I wasn't aware of it.

Thanks again for your help, and will appreciate any further
information on this last point of clarification regarding where R
references pixels (upper left or middle of pixel)?


Cheers, Lyndon

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

--
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