[gdal-dev] Writing a Tiff from R (was Writing a Tiff)

[Michael Sumner]

Hi,

This should work better. I'm not greatly familiar with the lower-level R 
functions for GDAL, but the following should be closer to what you wanted.

This assumes that "bm" is defined in the top level workspace, and the 
function can access it via R's scoping rules. You could use a similar 
approach to read the source in small pieces as well if need be. Please 
consult the rgdal documentation. I believe further discussion belongs on 
R-Sig-Geo as your problem likely has nothing to do with GDAL.

You may need to work with a column-wise version of indexing from your 
source matrix, depending on whether "bm" is aligned according to the 
image() conventions.

Cheers, Mike.

library(rgdal)
data(volcano)
bm <- volcano

##bm <- matrix(rnorm(1e6), nrow = 1000)  ## larger example

bm.writeTiff <- function(filename) {
       
        ## search for "bm" matrix in workspace (simplistic way to avoid 
copying)
        if (!(exists("bm") | is.matrix(bm)) ) stop("bm matrix not found ")
       
        # GTiff driver
        driver <- new('GDALDriver', 'GTiff')
       
       
        tds <- new("GDALTransientDataset",driver,nrow(bm),ncol(bm), 
type="Int32")

        for (row in 1:nrow(bm)) {
            ## note that leaving one index empty equates to specifying 
"all indices"
                putRasterData(tds, bm[row, ], offset= c(row - 1, 0)) ## 
R is 1-based, GDAL rows are 0-based
        }
       
        saveDataset(tds,filename)
       
        GDAL.close(tds)
}

bm.writeTiff("test.tif")


Henning Bredel wrote:
> On Mon, 2009-01-19 at 20:57 +1100, Michael Sumner wrote:
>   
>> Hello,
>>
>>     
>
> Hey Mike,
>
> thanks for your exhaustive answere :-)
>
>   
>> Note this recent R-Sig-Geo post about the "in-development" nature of the 
>> latest builds of rgdal - you need to let us know what versions you are 
>> using. Also, view the definition of writeGDAL() to see how the lower 
>> level functions are used.
>>
>> https://stat.ethz.ch/pipermail/r-sig-geo/2009-January/004830.html
>>
>>     
>
> I have read that message, but my version is not affected be the one
> Roger refers to:
>
> ,----------------------------------
>   
>> library(rgdal)
>>     
> Loading required package: sp
> Geospatial Data Abstraction Library extensions to R successfully loaded
> Loaded GDAL runtime: GDAL 1.5.2, released 2008/05/29
> GDAL_DATA: /usr/local/lib/R/site-library/rgdal/gdal
> Loaded PROJ.4 runtime: Rel. 4.6.0, 21 Dec 2007
> PROJ_LIB: /usr/local/lib/R/site-library/rgdal/proj
>   
>> sessionInfo()
>>     
> R version 2.7.1 (2008-06-23) 
> i486-pc-linux-gnu 
>
> locale:
> LC_CTYPE=en_US.UTF-8;LC_NUMERIC=C;LC_TIME=en_US.UTF-8;\
> LC_COLLATE=en_US.UTF-8;LC_MONETARY=C;LC_MESSAGES=en_US.UTF-8;\
> LC_PAPER=en_US.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;\
> LC_MEASUREMENT=en_US.UTF-8;LC_IDENTIFICATION=C
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods
> base     
>
> other attached packages:
> [1] rgdal_0.5-33 sp_0.9-29   
>
> loaded via a namespace (and not attached):
> [1] grid_2.7.1      lattice_0.17-20
> `----------------------------------
>
>   
>> BTW, you can simplify the writing of GTiffs with writeGDAL() -unless you 
>> are intending to develop your function to write to the file in smaller 
>> pieces. The code to do that can be as simple as this:
>>     
>
>
>   
>> library(rgdal)
>> ## read in matrix to bm
>> ## . . .
>> ## generate X/Y coordinates - assuming your matrix is aligned as image() 
>> expects (replace 0 and 1 as needed)
>> xx <- seq(0, length = nrow(bm), by = 1)
>> yy <- seq(0, length = ncol(bm), by = 1)
>> ##
>> writeGDAL(image2Grid(list(x = xx, y = yy, z = bm))
>>
>>     
>
> That snipped brought my machine »to its knees« .. but that is
> a problem of R itself (actually, I'm currently exploring).
>
>   
>> But, your R function works for me for the following simple test. As 
>> Frank mentioned, it probably comes down to a version problem or a file 
>> size limit (although that seems unlikely since your input matrix "bm" 
>> must be completely loaded in memory).
>>
>>  library(rgdal)
>> #Loading required package: sp
>> #Geospatial Data Abstraction Library extensions to R successfully loaded
>> #Loaded GDAL runtime: GDAL 1.6.0, released 2008/12/04
>> #Path to GDAL shared files: C:\inst\R\library/rgdal/gdal
>> #Loaded PROJ.4 runtime: Rel. 4.6.1, 21 August 2008
>> #Path to PROJ.4 shared files: C:\inst\R\library/rgdal/proj
>>
>> data(volcano)  # matrix 87x61
>> bm.writeTiff(volcano, "volc.tif")
>> #         used (Mb) gc trigger (Mb) max used (Mb)
>> #Ncells 305831  8.2     597831   16   407500 10.9
>> #Vcells 149194  1.2     786432    6   467409  3.6
>> tst <- readGDAL("volc.tif")
>>
>> #volc.tif has GDAL driver ltalk X6
>> #and has 87 rows and 61 columns
>>
>>     
>
> My image got dimension 8950 9790    7, could that dimension
> be the problem?
>
> Regards
>
>   Henning
>
>