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 Robert On Mon, Nov 15, 2010 at 3:51 AM, Roger Bivand <roger.biv...@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@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.biv...@nhh.no > > _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
