[R-sig-Geo] raster or rgdal: problem with geometry of tif + tifw file
Agustin Lobo
alobolistas at gmail.com
Mon May 9 14:51:54 CEST 2011
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