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.biv...@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.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 >>>> >>>> >>> >> >> > > -- > 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
