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.hijm...@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) les...@princeton.edu _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
