[R-sig-Geo] Pixel reference location?
Robert J. Hijmans
r.hijmans at gmail.com
Sun Nov 14 01:48:20 CET 2010
Hi Lyndon,
That clarifies what happens with 'crop' I think; but the different
georeferencing info you get from Arc and ENVI is puzzling. Your
question about georeferencing was clear, but my answer was not.
Another try: raster uses the extreme of the extreme cells for
georeferencing as in:
xoooo
ooooo
ooooo
i.e. (1,1), not (1.5, 1.5)
However, you say that Arc & ENVI put the coordinates half a cell
beyond where raster puts them (and now I understand where your
question came from), this is pretty serious as it suggest that they
are reading something different in the file. Could you send me the
file?
The cropping part:
to.crop <- raster(nrow=9872, ncol=6976, xmn = -283139.2, xmx =
1332896, ymn = -1172417, ymx = 1114495)
crop.img <- raster(nrow=6116, ncol=4200, xmn = -120863.9, xmx =
852092.8, ymn = -497022.6, ymx = 919787.6)
# What happens here:
cropped <- crop(to.crop, crop.img)
# is:
#1) get the intersection of to.crop and crop.img
e = intersectExtent(to.crop, crop.img)
#2) align the intersection with the raster that needs to be cropped
ee = alignExtent(e, to.crop)
cropped <- crop(to.crop, ee)
I suggested to use the extent function but I am now less sure about
that because there seems to be something else going on.
Best, Robert
On Sat, Nov 13, 2010 at 11:58 AM, Lyndon Estes <lestes at princeton.edu> wrote:
> 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).
It was just confusing to debug, if it all works it might be OK
>
> 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