A few months back I was getting some help with adapting R for doing tiled
processing of remote sensing data, and I made some serious headway, but
there is a lingering problem I keep scratching my head at. First off, a
quick recap of the base algorithm:
1) Let i = input raster, f() some function to apply to i, and o the output
raster, such that o = f(i)
2) i is read in "tiles", each tile read from the input we'll call i_tile (to
make it simplest, assume I'm just reading i one line at a time) -- this
process is performed using readGDAL(...offset=c(j, 0),...) where j is the
current line I'm reading. o_tile = f(i_tile), which is written using
*writeBin* (not writeGDAL).
3) This produces a flat binary file of the sort that Arc binary rasters or
ENVI rasters use.
So I produce the binary file, but the portion I'm getting a bit hung up on
is the production of the little text header which completes the file
(.hdr/.prj for an ESRI raster, .hdr for an ENVI raster, note that the .hdr
is different between arc and envi). For this example, let's try to make an
ENVI header. Remember i_tile is one of the readGDAL objects of a tiled
subset of image of filename i:
***
gp <- gridparameters(i_tile)
cellsize <- gp$cellsize
offset <- gp$cellcentre.offset
dims <- gp$cells.dim
dims[[2]]=rows
nbands <- 1
tempopen <- GDAL.open(i)
file_projection <- getProjectionRef(tempopen)
GDAL.close(tempopen)
# This was from Roger Bivand
d.drv <- new("GDALDriver", "ENVI")
tds.out <- new("GDALTransientDataset", driver = d.drv, rows = dims[2],cols
= dims[1], bands = nbands, type = "Float32")
gt <- c(offset[1] - 0.5 * cellsize[1], cellsize[1], 0.0,offset[2] +
(dims[2] -0.5) * cellsize[2], 0.0, -cellsize[2])
.Call("RGDAL_SetGeoTransform", tds.out, gt, PACKAGE = "rgdal")
.Call("RGDAL_SetProject", tds.out, file_projection, PACKAGE = "rgdal")
fn <- tempfile()
saveDataset(tds.out, fn)
GDAL.close(tds.out)
GDAL.close(d.drv)
temporary_dir=dirname(fn)
outdir=dirname(outfile)
fnhdr<-as.character(paste(basename(fn),outfile_hdr,sep=''))
outhdr<-as.character(paste(basename(rasteroutput),outfile_hdr,sep=''))
setwd(temporary_dir)
file.rename(fnhdr,outhdr)
file.copy(outhdr,outdir,overwrite=TRUE)
***
So this WORKS the first time I run it, but the SECOND time it is run
(without quitting out of R) I get the following error:
d.drv <- new("GDALDriver", "ENVI")
Error in .local(.Object, ...) : No driver registered with name: ENVI
It seems like there is some lingering open driver or something someplace,
but for the life of me I can't tell where it is. You can see that I did a
GDAL.close(d.drv) but that doesn't seem to be sufficient. Thoughts?
By the way, the full algorithm would deal with RS data of an arbitrary
number of samples, lines and bands (e.g. hyperspectral processing would be
completely feasible). I'm happy to send the script to the list once the
bugs are worked out! Thanks!
--j
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
