raster or rgdal: problem with geometry of tif + tifw file
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:
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:
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()),
writeRaster(round(testcirpred),file="/media/TRANSCEND2/VegcamPro/MATA20090729/Proc2/delme.tif", format="GTiff", overwrite=TRUE,datatype='INT2U')
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>:
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:
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:
test= brick("/media/TRANSCEND2/VegcamPro/MATA20090729/TIFFOri/SDIM0114.tif")
show(test)
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:
writeRaster(test,file="test.tif", format="GTiff", overwrite=TRUE,datatype='FLT8S')
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
sessionInfo()
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
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo