On Thu, Jun 11, 2009 at 3:50 PM, Roger Bivand <Roger.Bivand@nhh.no> wrote:

> On Thu, 11 Jun 2009, Rainer M Krug wrote:
>
>  Hi
>>
>> I spend the whole day with it, trying grass, gdal, R, but was not
>> successful
>> (and I know it can be done...).
>>
>> I have a shape file and a raster file, both in the same projection and in
>> thew same area (shape file completely overlapping the raster) and I want
>> to
>> "burn" the values from one column in the attribute table of the shape file
>> (DNSTY) into the raster.
>>
>> I am doing the following:
>>
>> ####
>> library(sp)
>> library(maptools)
>> library(raster) # <- this might work as well. I also have the shape file
>> as
>> a .csv file, including coordinates and DNSTY
>>
>> s <- readShapePoints("./plants.shp")
>> r <- readAsciiGrid("./plants.asc")
>>
>> ## and now?
>> ####
>>
>> I seem to remember that it works by using overlay(), but I don't remember
>> how.\
>>
>
> Please see the thread started by Karla Falk yesterday. overlay() will tell
> you which raster cells the points fall in, but you cannot just "burn" across
> unless you know definitely that hit raster cells have only one hit,
> otherwise you'll need tapply (or table for counting) to tally the point
> values by hit cell. I think that is the step you need to add. It'll take
> some attempts to find a recipe that works for your specific case, it usually
> does, but then works.
>

Thanks for your tip - I found a solution. Below, I'll outline it:

library(maptools)
r <- readAsciiGrid("MASK.asc")
r[[1]] <- NA
p <- readShapePoints("plants.shp")

## check if more then one point falls in a single cell
to <- table(rownames(o@data))
sum(to > 1)
0 ## OK - not the case, so we can continue.

## set the values of the cells in the raster r equal to the values of the
points in p
## points are identified by the rowname of the data
r[[1]][rownames(r@data) %in% rownames(o@data)] <- p$DNSTY

writeAsciiGrid(r, "plants.asc")

and that's it.

To make life easier, I think a function in sp, which takes a
SpatialPointsDataFrame p, a SpatialGridDataFrame g and a function FUN used
to aggregate the values when more then one point is in a single cell, would
be a good idea.

Thanks again for your help,

Rainer


> Roger
> SpatialGridDataFrame
>
>  Also, it should be working with pointsToRaster() from the package raster,
>> but I also was not successful.
>>
>> Cheers,
>>
>> Rainer
>>
>>
>>
> --
> Roger Bivand
> Economic Geography Section, Department of Economics, Norwegian School of
> Economics and Business Administration, Helleveien 30, N-5045 Bergen,
> Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
> e-mail: Roger.Bivand@nhh.no
>
>


-- 
Rainer M. Krug, Centre of Excellence for Invasion Biology, Stellenbosch
University, South Africa

	[[alternative HTML version deleted]]

