[R-sig-Geo] raster or rgdal: problem with geometry of tif + tifw file
Robert J. Hijmans
r.hijmans at gmail.com
Wed May 11 21:11:17 CEST 2011
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
>>>
>>
>
More information about the R-sig-Geo
mailing list