[R-sig-Geo] Georerferncing in Arc vs rgdal (and raster): AREA_OR_POINT=Point
Roger Bivand
Roger.Bivand at nhh.no
Mon Nov 15 10:51:56 CET 2010
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 at 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.Bivand at nhh.no
More information about the R-sig-Geo
mailing list