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

Roger Bivand Roger.Bivand at nhh.no
Tue Nov 16 09:26:58 CET 2010


On Mon, 15 Nov 2010, Roger Bivand wrote:

> On Mon, 15 Nov 2010, Robert J. Hijmans wrote:
>
>> 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
>
> Good, thanks Robert.
>
> There is already an indication that someone from ESRI has answered, not on 
> your ticket, but on the response you cite. I've replied to 
> http://trac.osgeo.org/gdal/ticket/3838, so we'll see what happens. Until it 
> is sorted out between ESRI and GDAL, I suggest that we wait - Franks 
> immediate response that this should (if necessary) be handled by 
> GDALDataset::GetGeoTransform() is sensible.

Update:

See http://trac.osgeo.org/gdal/wiki/rfc33_gtiff_pixelispoint for current 
status. Briefly, the legacy GeoTIFF documentation said one thing and the 
GDAL documentation and implementation said and did another. GDAL behaviour 
will change from 1.6.4, 1.7.4 and 1.8.0 if AREA_OR_POINT has been set to 
its non-default value of POINT, and the new GTIFF_POINT_GEO_IGNORE 
configuration option is TRUE (default for 1.6 and 1.7).

ESRI were following the GeoTIFF documentation, GDAL wasn't, and the RFC 
explains:

"Traditionally GDAL has treated this flag as having no relevance to the 
georeferencing of the image despite disputes from a variety of other 
software developers and data producers. This was based on the authors 
interpretation of something said once by the GeoTIFF author".

Roger


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

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