[R-sig-Geo] Georerferncing in Arc vs rgdal (and raster): AREA_OR_POINT=Point
Robert J. Hijmans
r.hijmans at gmail.com
Tue Nov 16 15:25:11 CET 2010
And I was wrong to think that Arc corrected incorrectly. I now think
they got it right. I have the same correction in raster 1.6-20; with a
warning, partly to remember to turn this correction off again when it
is fixed in gdal. Best, Robert
On Tue, Nov 16, 2010 at 2:26 AM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> 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