[R-sig-Geo] Fwd: writeRaster problems

[Robert J. Hijmans]

Dear Seth,

I think the problem is that you try to overwrite values in an existing
file. That is not directly possible (with raster anyway). You have to
start a new fil. Below are some examples. Also see
vignette('functions', 'raster')

#1 if you can read the whole raster into memory
p1a <- rep(0, ncol(rout))
rin <- raster(system.file("external/test.grd", package="raster"), values=TRUE)
rin[1:ncol(rin)] = p1a
rout <- writeRaster(rin, filename="write1.rst", format="IDRISI",overwrite=TRUE)

#2 row by row (and slow)
rin <- raster(system.file("external/test.grd", package="raster"))
rout <- raster(rin)
p1a <- rep(0, ncol(rout))
rout <- setValues(rout, p1a, 1)
rout <- writeRaster(rout, filename="write2.rst", format="IDRISI",overwrite=TRUE)
for (r in 2:nrow(rout)) {
       rin <- readRow(rin, r)
       rout <-setValues(rout, values(rin), r)
       rout <- writeRaster(rout, filename="write5.rst",
format="IDRISI",overwrite=TRUE)
}


# 3 block by block (much faster); using the rgdal driver
(format='RST') so that you
# first copy all values and then overwrite rows at any location
rin <- raster(system.file("external/test.grd", package="raster"))
rout <- raster(rin)
bs <- blockSize(rout)
rout <- writeStart(rout, filename="write3.rst", format="RST", overwrite=TRUE)
for (i in 1:bs$n) {
       v = getValuesBlock(rin, row=bs$row[i], nrows=bs$size)
       writeValues(rout, v, bs$row[i])
}
p1a <- rep(0, ncol(rout))
writeValues(rout, p1a, 1) # replace a row
rout <- writeStop(rout)


Robert



On Wed, May 5, 2010 at 6:52 PM, Seth J Myers <sjmyers at syr.edu> wrote:
> Hi,
>
> I'm running R 2.11.0 on a 32 bit Windows Vista machine.  I have new versions of all relevant libraries installed.  I have been reading from IDRISI .rst files and this works fine using raster() and getValues().  I have passed the values from the .rst files through another R object and would like to write them as an .rst file.  In the following code, I am attempting to set the 1st row values of a RasterLayer object and then write these to a .rst file.  All raster files I am working with have 3,700 columns and 4,203 rows and are Connecticut StatePlane coordinate system which is equivalent to lambert's conformal conic proj.
>
> #extract 2nd column of a matrix which are the values of interest for the 1st row
> p1a<-p1[c(1:3700),c(2)]
> p1am<-as.matrix(p1a)
> #create RasterLayer object from a .rst file that was created in IDRISI beforehand with all 0s (real nums) for cell values
> put<-raster("C:\\rforest\\data\\rst85\\write2.rst",values=FALSE)
> #set 1st row to equal the values of interest
> put<-setValues(put,p1am,1)
> #write to a .rst file that I created in IDRISI beforehand with all 0s (real nums)
> writeRaster(put, filename="C:\\rforest\\data\\rst85\\write5.rst", format="IDRISI",overwrite=TRUE)
>
>
> PROBLEM:
>
> Until the writeRaster function call, everything is working.  I can use getValues to look at the 1st row of "put", and this equals p1am values.  I can look at the 2nd row of "put" and they are all zeros.  The writeRaster function does not give an error but, when opened in IDRISI, the legend seems correct but the values are expressed in scientific notation and are MUCH too small.  Also, the 1st row values are copied for all rows, when I expect all but row=1 to be equal to 0.  I tried, prior to opening in IDRISI, to create a RasterLayer object from the output raster from writeRaster above, to see if R gives the result wanted.  This results in the following error message and report on the RasterLayer object created.
>
>> put2<-raster("C:\\rforest\\data\\rst85\\write5.rst",values=TRUE)
>> put3<-getValues(put2,1)
> Error:
>        GDAL Error 3: Can't read(C:\rforest\data\rst85\write5.rst) block with X offset 0 and Y offset 0.
> No error
> Error in values(readRows(x, row, nrows)) :
>  error in evaluating the argument 'x' in selecting a method for function 'values'
>> put2
> class       : RasterLayer
> filename    : C:\rforest\data\rst85\write5.rst
> nrow        : 4203
> ncol        : 3700
> ncell       : 15551100
> min value   : 0
> max value   : 0
> projection  : +proj=lcc +lat_1=41.86666666666667 +lat_2=41.2 +lat_0=40.83333333333334 +lon_0=-72.75 +x_0=304800.6096 +y_0=152400.3048 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0
> xmin        : 221895.5
> xmax        : 334671.7
> ymin        : 164342.4
> ymax        : 292450.1
> xres        : 30.48006
> yres        : 30.48006
>
>
> I have tried this entire process starting NOT with blank 0 rasters to create my SpatialLayer object (a shortcut I thought may be a problem) but instead creating rasters by specifying projections etc on non-used filenames.  The results are the same.
>
> Thanks,
> Seth
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>