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]
}
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
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.biv...@nhh.no
