Skip to content

raster or rgdal: problem with geometry of tif + tifw file

5 messages · Agustin Lobo, Robert J. Hijmans

Original
Agustin Lobo
# 
Hi!

I read a tif fle (with a companion tifw file) for which gdalinfo states:

alobo at alobo-laptop:/media/TRANSCEND2/VegcamPro/MATA20090729/Proc2$
gdalinfo SDIM0114.tif
Warning 1: TIFFReadDirectory:Unknown field with tag 37717 (0x9355) encountered
Driver: GTiff/GeoTIFF
Files: SDIM0114.tif
Size is 2640, 1760
Coordinate System is `'
GeoTransform =
  445548.7745816645, -0.109570548549058, 0.216912817625298
  4628418.279379116, 0.216912817625298, 0.109570548549058
Metadata:
  TIFFTAG_SOFTWARE=SIGMA Photo Pro 3.5.2.0000
  TIFFTAG_DATETIME=2010:02:27 15:40:55
  TIFFTAG_XRESOLUTION=180
  TIFFTAG_YRESOLUTION=180
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  445548.775, 4628418.279)
Lower Left  (  445930.541, 4628611.124)
Upper Right (  445259.508, 4628990.929)
Lower Right (  445641.275, 4629183.773)
Center      (  445595.025, 4628801.026)


Then I read it in R:
class       : RasterBrick
dimensions  : 1760, 2640, 3  (nrow, ncol, nlayers)
resolution  : 0.1095705, 0.1095705  (x, y)
extent      : 445548.8, 445838, 4628418, 4628611  (xmin, xmax, ymin, ymax)
projection  : NA
values      : /media/TRANSCEND2/VegcamPro/MATA20090729/TIFFOri/SDIM0114.tif
min values  : 0 0 0
max values  : 65535 65535 65535

and write it back to disk:
But gdalinfo states a different bounding box for the new file:
alobo at alobo-laptop:/media/TRANSCEND2/VegcamPro/MATA20090729/Proc2$
gdalinfo test.tif
Driver: GTiff/GeoTIFF
Files: test.tif
Size is 2640, 1760
Coordinate System is `'
Origin = (445548.774581664009020,4628611.233115110546350)
Pixel Size = (0.109570548549241,-0.109570548548469)
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  445548.775, 4628611.233)
Lower Left  (  445548.775, 4628418.389)
Upper Right (  445838.041, 4628611.233)
Lower Right (  445838.041, 4628418.389)
Center      (  445693.408, 4628514.811)
.../...

Somehow the geometric properties are modified in R

I've made a single geotif file with gdal_transform (no tifw created)
and read it in R in the same way and get the same result, so the
problem does not come
from using a tifw file as I initially thought.

Perhaps a bug in brick() ?

Thanks

Agus

Files:
http://dl.dropbox.com/u/3180464/SDIM0114.tif
http://dl.dropbox.com/u/3180464/SDIM0114.tifw
R version 2.13.0 (2011-04-13)
Platform: i486-pc-linux-gnu (32-bit)

locale:
 [1] LC_CTYPE=en_US.UTF-8          LC_NUMERIC=C
LC_TIME=en_US.UTF-8
 [4] LC_COLLATE=en_US.UTF-8        LC_MONETARY=en_US.UTF-8
LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=en_US.UTF-8          LC_NAME=en_US.UTF-8
LC_ADDRESS=en_US.UTF-8
[10] LC_TELEPHONE=en_US.UTF-8      LC_MEASUREMENT=en_US.UTF-8
LC_IDENTIFICATION=en_US.UTF-8

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets
methods   base

other attached packages:
[1] gridExtra_0.7 ggplot2_0.8.8 rgdal_0.6-31  raster_1.8-15 sp_0.9-72
   reshape_0.8.3
[7] plyr_1.2.1    proto_0.3-8   rkward_0.5.6

loaded via a namespace (and not attached):
[1] lattice_0.19-26 tools_2.13.0
Robert J. Hijmans
# 
Hi Agus,

The problem is that you raster is rotated (GeoTransform[3] and [5] are not zero)

GeoTransform =
 445548.7745816645, -0.109570548549058, 0.216912817625298
 4628418.279379116, 0.216912817625298, 0.109570548549058

raster ignored the rotation. I have fixed that (in version 1-8-19) by
throwing an error for such a file.

readGDAL will read the file, but as points (SpatialPointsDataFrame),
which you could interpolate to a non-rotated raster. I do not plan to
support such files in raster, but I might include a function for
memory-safe transformation to a non-rotated file.

Thanks for catching & reporting this,

Robert
On Sun, May 8, 2011 at 9:50 AM, Agustin Lobo <alobolistas at gmail.com> wrote:
Agustin Lobo
# 
Robert,

I think a warning is better than an error. There are many operations
that can be performed
ignoring the rotation, it should be up to the user.

In my case the problem is actually created at writing (writeRaster()),
because the new tif file gets wrong  geometric information that has
priority over the one provided by the tfw file.

I circumvent the problem by using
listgeo -no_norm SDIM0114.tif > original.geo
in a directory where the tfw file is not present. Then I apply this
geometric information to the output from R:
geotifcp -g original.geo delme.tif SDIM0114L1.tif
cp SDIM0114.tfw SDIM0114L1.tfw

more or less as suggested in
http://remotesensing.org/geotiff/faq.html

It would be good being able to just omit the geometric info for the
output tif file in writeRaster(), so the original tfw file would be
directly valid for the output file
to be displayed in QGIS, thus avoiding the need of the fix with
listgeo and geotifcp

Agus

2011/5/8 Robert J. Hijmans <r.hijmans at gmail.com>:
Robert J. Hijmans
# 
Hi Agus,

Thanks for the suggestions. I'll make it a (stern) warning. I am also
going to store the rotation info add build in some support for such
files.

I believe you can write a "plain tiff" file by using

options="PROFILE=BASELINE"

as in:

x = writeRaster(r, 'file.tif', options="PROFILE=BASELINE")

Robert
On Mon, May 9, 2011 at 5:51 AM, Agustin Lobo <alobolistas at gmail.com> wrote:
2 days later
Robert J. Hijmans
# 
Hi Agus,

I have added  some support for rotated rasters to the raster package
(>= 1.8-20).

- There is a warning
- Most functions that assume rectangularity, will now fail, with a
clear warning.
- Some of these functions, such as plot, do not need to fail, but need
to be fixed with a few more lines of code. I will do that, little by
little, but let me know if you find an important one.
- If you write results to a new geotiff the rotation parameters should
now be preserved correctly
- A new function "rectify" to get rid of the rotation

library(raster)
beginCluster() # speeds things up tremendously, if you have a multi-core machine
b <- brick('SDIM0114.tif')
br <- rectify(b, progress='window')
plotRGB(br)

Robert
On Mon, May 9, 2011 at 11:17 AM, Robert J. Hijmans <r.hijmans at gmail.com> wrote: