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

Mon Nov 15 16:11:12 CET 2010

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.

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