[R-sig-Geo] writeRaster problems
Robert J. Hijmans
r.hijmans at gmail.com
Thu May 6 20:11:11 CEST 2010
> 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.pebesma at 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.hijmans at 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 <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
>>>>
>>>
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at 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.pebesma at wwu.de
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
More information about the R-sig-Geo
mailing list