[R-sig-Geo] raster or rgdal: problem with geometry of tif + tifw file
Robert J. Hijmans
r.hijmans at gmail.com
Mon May 9 20:17:15 CEST 2011
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