> Please note that in order to get the usual map view for this data set,
> such that 1 km NS equals 1 km EW, you need to add:
>
> plot(r, asp=1)
And this is now the default behavior ( raster >= 1.0.8 ). Thanks for
pointing this out (again).
On Wed, May 5, 2010 at 11:18 PM, Edzer Pebesma
<edzer.pebe...@uni-muenster.de> wrote:
> On 05/06/2010 07:26 AM, Robert J. Hijmans wrote:
>> 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)
>
> Please note that in order to get the usual map view for this data set,
> such that 1 km NS equals 1 km EW, you need to add:
>
> plot(r, asp=1)
>
>> 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
>
> --
> Edzer Pebesma
> Institute for Geoinformatics (ifgi), University of Münster
> Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
> 8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de
> http://www.52north.org/geostatistics e.pebe...@wwu.de
>
> _______________________________________________
> 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
