Georerferncing in Arc vs rgdal (and raster): AREA_OR_POINT=Point
Hi Robert, In case this adds a piece to the puzzle, here is the version information : R version 2.11.1 (2010-05-31) x86_64-apple-darwin9.8.0 raster_1.5-16 sp_0.9-71 Thanks, Lyndon
On Sun, Nov 14, 2010 at 2:27 AM, Robert J. Hijmans <r.hijmans at gmail.com> wrote:
List, Lyndon Estes asked questions about georeferencing and the use of 'crop' in raster, while pointing out differences in georeferencing between Arc (and ENVI) vs raster (and rgdal). I start a new thread to focus on the georeferencing issue. I think there is something going seriously wrong. GDAL supports metadata ( http://www.gdal.org/gdal_datamodel.html ). One of the items is "AREA_OR_POINT" which may be either "Area" (the default) or "Point". This "Indicates whether a pixel value should be assumed to represent a sampling over the region of the pixel or a point sample at the center of the pixel. This is not intended to influence interpretation of georeferencing which remains area oriented." The file Lyndon send me ("MOD13Q1.A2005225.h20v11.mosaic.NDVI.tif") has AREA_OR_POINT=Point (reported by rgdal, see further below). This is ignored by rgdal and raster (as it should). However, it appears that ArcMap version 9.3 (and ENVI?) does not ignore this flag and changes the georeference.
#this is what arc reports arc <- extent(-283255.039878, 1332779.71537, -1172300.9282, 1114610.64058) extent(arc)
class ? ? ? : Extent xmin ? ? ? ?: -283255.0 xmax ? ? ? ?: 1332780 ymin ? ? ? ?: -1172301 ymax ? ? ? ?: 1114611
# make a RasterLayer of the file
r1 <- raster('MOD13Q1.A2005225.h20v11.mosaic.NDVI.tif')
extent(r1)
class ? ? ? : Extent xmin ? ? ? ?: -283139.2 xmax ? ? ? ?: 1332896 ymin ? ? ? ?: -1172417 ymax ? ? ? ?: 1114495
# different by half a cell res(r1)
[1] 231.6564 231.6564
The weird thing is that Arc seems to correct (which I think it should not) and then makes a mistake in the correction. This is how you can get the georeference that Arc has:
xmin(r1) - 0.5 * xres(r1)
[1] -283255.0
xmax(r1) - 0.5 * xres(r1)
[1] 1332780
ymin(r1) + 0.5 * yres(r1)
[1] -1172301
ymax(r1) + 0.5 * yres(r1)
[1] 1114611
I.e., shifting the xmin AND xmax half a cell to the left, and the ymin AND ymax half a cell up. It makes no sense to do that as you would want to shift the xmin to the left and the xmax to the right to go from center-of-cell georeferencing to extreme-of -cell georeferencing. With recent versions of raster, you can do this correction like this ( not documented, sorry ):
r2 <- raster('MOD13Q1.A2005225.h20v11.mosaic.NDVI.tif', fixGeoref=TRUE)
extent(r2)
class ? ? ? : Extent xmin ? ? ? ?: -283255.1 xmax ? ? ? ?: 1333011 ymin ? ? ? ?: -1172533 ymax ? ? ? ?: 1114611
And now we have the same xmin and ymax as Arc, but the xmax and ymin are different (and I believe raster is right here) When I save the data as a new geotiff file, and by default, with no AREA_OR_POINT=Area r3 <- writeRaster(r1, filename='test.tif') Arc and rgdal/raster agree again. That suggests that it is this flag that matters, but perhaps there is something else in the file to which Arc responds? Does anyone have any insights, thoughts? Did I miss something simple? Or does Arc have a "standard" (perhaps shared with many others?) that we should be aware of? More details below Robert
GDALinfo('MOD13Q1.A2005225.h20v11.mosaic.NDVI.tif')
rows ? ? ? ?9872 columns ? ? 6976 bands ? ? ? 1 origin.x ? ? ? ?-283139.2 origin.y ? ? ? ?-1172417 res.x ? ? ? 231.6564 res.y ? ? ? 231.6564 ysign ? ? ? -1 oblique.x ? 0 oblique.y ? 0 driver ? ? ?GTiff 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 file ? ? ? ?MOD13Q1.A2005225.h20v11.mosaic.NDVI.tif apparent band summary: ?GDType ? Bmin ?Bmax Bmean Bsd 1 ?Int16 -32768 32767 ? ? 0 ? 0 Metadata: TIFFTAG_SOFTWARE=HEG-Modis Reprojection Tool ?Nov 4, 2004 AREA_OR_POINT=Point Warning message: statistics not supported by this driver
# this is how the I fix the georeferences to go from centre to extreme
based, with x being a Raster object
? ? ? ?if (fixGeoref) {
? ? ? ? ? ? ? ?xx <- x
? ? ? ? ? ? ? ?nrow(xx) <- nrow(xx) - 1
? ? ? ? ? ? ? ?ncol(xx) <- ncol(xx) - 1
? ? ? ? ? ? ? ?rs <- res(xx)
? ? ? ? ? ? ? ?xmin(x) <- xmin(x) - 0.5 * rs[1]
? ? ? ? ? ? ? ?xmax(x) <- xmax(x) + 0.5 * rs[1]
? ? ? ? ? ? ? ?ymin(x) <- ymin(x) - 0.5 * rs[2]
? ? ? ? ? ? ? ?ymax(x) <- ymax(x) + 0.5 * rs[2]
? ? ? ?}
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