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
_______________________________________________
gdal-dev mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
