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

[Lyndon Estes]

Hi Robert,

In case this adds a piece to the puzzle, here is the version information :

R version 2.11.1 (2010-05-31)
x86_64-apple-darwin9.8.0

raster_1.5-16
sp_0.9-71

Thanks, Lyndon


On Sun, Nov 14, 2010 at 2:27 AM, Robert J. Hijmans <r.hijmans at gmail.com> 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."
>
> 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]
>        }
>



-- 
Lyndon Estes
Research Associate
Woodrow Wilson School
Princeton University
+1-609-258-2392 (o)
+1-609-258-6082 (f)
+1-202-431-0496 (m)
lestes at princeton.edu