raster package: overlay changes geometry
Hi Agus,
This is probably caused by unexpected values in your input TIF file.
The Upper Left Y is lower then the Lower Left Y...
Upper Left ( 0.0, 0.0)
Lower Left ( 0.0, 1760.0)
Upper Right ( 2640.0, 0.0)
Lower Right ( 2640.0, 1760.0)
I do not know if this is "legal" (I sure hope not); but it would be
useful to have your input file (if you can send it off-list) to see
how to best catch this (and either deal with it correctly or give an
error).
I cannot reproduce this problem with other files
Here is an example that works fine:
library(raster)
fn1 <- system.file("pictures/Rlogo.jpg", package="rgdal")
b <- brick(fn1)
b <- readAll(b)
# create a three band TIF file as in your example
b2 <- writeRaster(b, filename='input.tif', filetype='GTiff', overwrite=TRUE)
fn2 <- 'input.tif'
r1 <- raster(fn2,band=1)
r2 <- raster(fn2,band=2)
r3 <- raster(fn2,band=3)
f <- function(x,y,z,coefs= t(matrix(1:3)))
coefs[1,1]*x+coefs[1,2]*y+coefs[1,3]*z
o <- overlay(r1,r2,r3,fun=f,filename="test.tif", filetype='GTiff',
overwrite=TRUE, progress='tcltk')
extent(b)
extent(b2)
extent(r1)
extent(o)
GDALinfo(fn1)
GDALinfo(fn2)
GDALinfo('test.tif')
# By the way, this would be more efficient
b <- brick(fn1)
overlay(b,fun=f,filename="test.tif", filetype='GTiff', overwrite=TRUE)
Robert
On Tue, Oct 27, 2009 at 3:44 AM, Agustin Lobo <alobolistas at gmail.com> wrote:
Hi!
I've done:
r1 <- raster("SDIM0160.tif",band=1)
r2 <- raster("SDIM0160.tif",band=2)
r3 <- raster("SDIM0160.tif",band=3)
f <- function(x,y,z,coefs=calcoefs) coefs[1,1]*x+coefs[1,2]*y+coefs[1,3]*z
overlay(r1,r2,r3,fun=f,filename="SDIM0160_miL1.tif"filetype='GTiff')
The problem is that the resulting raster does not have the same geometry as
the input one.
Input:
$ gdalinfo SDIM0160.tif
Warning 1: TIFFReadDirectory:Unknown field with tag 37717 (0x9355)
encountered
Driver: GTiff/GeoTIFF
Files: SDIM0160.tif
Size is 2640, 1760
Coordinate System is `'
Metadata:
?TIFFTAG_SOFTWARE=SIGMA PhotoPro 2.4.0.3198
?TIFFTAG_DATETIME=2009:09:22 15:56:48
?TIFFTAG_XRESOLUTION=180
?TIFFTAG_YRESOLUTION=180
?TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
Image Structure Metadata:
?INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ?( ? ?0.0, ? ?0.0)
Lower Left ?( ? ?0.0, 1760.0)
Upper Right ( 2640.0, ? ?0.0)
Lower Right ( 2640.0, 1760.0)
Center ? ? ?( 1320.0, ?880.0)
Band 1 Block=2640x1 Type=Byte, ColorInterp=Red
Band 2 Block=2640x1 Type=Byte, ColorInterp=Green
Band 3 Block=2640x1 Type=Byte, ColorInterp=Blue
Output:
$ gdalinfo SDIM0160_miL1.tif
Driver: GTiff/GeoTIFF
Files: SDIM0160_miL1.tif
Size is 2640, 1760
Coordinate System is `'
Origin = (0.000000000000000,1761.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Image Structure Metadata:
?INTERLEAVE=BAND
Corner Coordinates:
Upper Left ?( ? ? ? 0.000, ? ?1761.000)
Lower Left ?( ? 0.0000000, ? 1.0000000)
Upper Right ( ? ?2640.000, ? ?1761.000)
Lower Right ( ? ?2640.000, ? ? ? 1.000)
Center ? ? ?( ? ?1320.000, ? ? 881.000)
Band 1 Block=2640x1 Type=Float32, ColorInterp=Gray
Could this be fixed?
Thanks,
Agus
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo