[R-sig-Geo] Georerferncing in Arc vs rgdal (and raster): AREA_OR_POINT=Point

[Robert J. Hijmans]

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