Seth,
I think I overlooked what is perhaps the easiest way to do this, using
a standard replacement approach:
library(raster)
r <- raster(system.file("external/test.grd", package="raster"))
p1a <- rep(100, ncol(r))
r[1:ncol(r)] <- p1a
r <- writeRaster(r, filename='test.rst', overwrite=TRUE)
plot(r)
Robert
On Wed, May 5, 2010 at 10:05 PM, Robert J. Hijmans <r.hijm...@gmail.com> wrote:
> 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 <sjmy...@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@stat.math.ethz.ch
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
