Georerferncing in Arc vs rgdal (and raster): AREA_OR_POINT=Point

Thanks Roger,

It would indeed be good to have a publicly available file for this.
But I now think that GDAL gets it wrong because the geotiff
specification ( see
http://www.remotesensing.org/geotiff/spec/geotiff2.5.html )
distinguishes between "PixelIsArea" and "PixelIsPoint" (Pixel is
point). The specification also says that these has implications for
georeferencing (shifting half a cell).

If AREA_OR_POINT=Point reflects "PixelIsPoint" then GDAL should not
ignore that (and if it does, we should make the correction ourselves).
Perhaps something (adjustment of coordinates) happens under the
gdal-hood but that does not seem to be the case.

While ArcMap corrects, it appears that they do not correct correctly,
by shifting the whole structure 0.5 cells to the left and up. It seems
to me that it should be that the upper left corner needs to shift 0.5
cells up and left, but the lower right corner needs to shift 0.5 cells
down and to the right.

I submitted a ticket to GDAL http://trac.osgeo.org/gdal/ticket/3837
and a response here: http://trac.osgeo.org/gdal/ticket/3838

Robert



On Mon, Nov 15, 2010 at 3:51 AM, Roger Bivand <Roger.Bivand at nhh.no> wrote:

On Sat, 13 Nov 2010, Robert J. Hijmans 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."

Exactly. The metadata records and reports the support of the grid
representation, that is whether the data "represent" only the central point
of the pixel (like a bore core), or "represent" the surface area of the
pixel (like a remote sensing pixel). In GTiff and GDAL more generally, it
has nothing whatsoever to do with GDALDataset::GetGeoTransform(), which
returns the upper left corner of the upper left pixel, and pixel width and
height for a dataset, and is where the georeferencing happens.

For example, users might like to set the metadata to "Area" for the output
of block kriging, and "point" for ordinary kriging, to record the support of
the data.

Whether ESRI adheres to this simple distinction is unknown (I have no such
software, even if I was motivated to check, which I am afraid isn't the
case). It is possible that they, unlike GDAL, have two raster cell models,
which do not fit the simple GDALDataset::GetGeoTransform() model, and are
using this support metadata item in a non-standard way to switch between
models. Software should not try to handle (too many cases of) odd behaviour
- these remain the user's responsibility to handle. If
GDALDataset::GetGeoTransform() from ESRI-generated (which version? only
9.3?) GTiff give non-standard results depending on the setting of this
metadata item, this isn't our problem, it is ESRI's problem.

Can someone please post a full set of use cases (files someone sent someone
are no use) with complete lineages (where the data started and what they've
been through along the way), preferably with copies of all the rasters all
the way along the workflow? Then maybe someone with access to a range of
versions of ESRI software can undertake to do due diligence. In addition, we
can ask the gdal-dev list to comment - although I think that the
documentation cited by Robert is conclusive - this has nothing to do with
georeference.

Roger

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

--
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no