[R-sig-Geo] Pixel reference location?
Lyndon Estes
lestes at princeton.edu
Sat Nov 13 20:58:23 CET 2010
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
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at stat.math.ethz.ch
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
--
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
More information about the R-sig-Geo
mailing list