[R-sig-Geo] Create raster layer that counts SpatialPoint occurrence within cells

[Mathieu Basille]

Hi Andrew,

The answer from Robert is very neat! Alternatively, you may also have a 
look at `count.points` in the package adehabitatMA (which was developed for 
animal locations!). Here is a working example from the vignette (p.17):

library(adehabitatMA)

## Load the data
data(lynxjura)

## `map` is a SpatialPixelsDataFrame
map <- lynxjura$map
class(map)

## `locs` is a SpatialPointsDataFrame
locs <- lynxjura$loc
class(locs)

## Plot the map + points
image(map, 1)
points(locs, pch=3)

## Count the number of points per pixel and map the result
cp <- count.points(locs, map)
image(cp)

## `cp` is a SpatialPixelsDataFrame
class(cp)

Hope this helps,
Mathieu.


Le 11/09/2014 19:01, Andrew Vitale a écrit :
> Hello all,
>
> I have a Python script that searches Twitter for tweets that meet certain
> criteria and records the coordinates of those tweets.  This leaves me with
> a CSV that can then easily be read into R as a SpatialPoints object.
>
> Once I have the SpatialPoints, I would like to make a raster layer that has
> the count of tweet locations that occur within each cell of that raster
> layer.
>
> I have a working script, but I feel like my way of doing this is very
> inefficient.  Does anyone know a better way to produce a map similar to
> this?  I'm also interested in other suggestions of how to visualize these
> data.
>
> The example below shows my approach to this problem.  I believe the weak
> point of the code is my for loop.  I would like to increase efficiency, as
> I plan to use this approach on much larger datasets.  Thanks for your
> advice!
>
> -Andrew Vitale
>
> The example:
>
>
> library(raster)
>
> ## Create a set of random long/lat coordinates
> ## that are centred around a "metro" area.
> ## These data mimic the CSV of tweet coords I read in to R
> n = 500
> set.seed(232)
> x = rnorm(n, -96, 5)
> y = rnorm(n, 32, 5)
> xy = cbind(x, y)
>
> ## Create a SpatialPoints object of the "tweets"
> tweets = SpatialPoints(xy)
> proj4string(tweets) = CRS('+proj=longlat +datum=WGS84')
>
> ## Create an empty template raster for rasterizing the points and
> ## create an empty list to store each rasterized point
> r = raster(nrow=90, ncol=180)
> l = list()
>
> ## Loop through the tweets one point (row) at a time and create
> ## a raster with the value of 1 at the location of the tweet.
> ## Store the raster layers created for each tweet in a list.
> ## This loop is the slow part of my code. ##
> for(i in 1:length(tweets)){
>    p = tweets[i, ]
>    l[[i]] = rasterize(p, r)
> }
>
> ## Create a raster stack from the list
> s = stack(l)
>
> ## Sum the stack layers together ignoring NA values, which
> ## produces a raster layer of counts.  Set 0 values to NA
> counts = sum(s, na.rm=TRUE)
> counts[counts == 0] = NA
>
> ## Plot the raster layer
> plot(counts)
>

-- 

~$ whoami
Mathieu Basille
http://ase-research.org/basille

~$ locate --details
University of Florida \\
Fort Lauderdale Research and Education Center
(+1) 954-577-6314

~$ fortune
« Le tout est de tout dire, et je manque de mots
Et je manque de temps, et je manque d'audace. »
  -- Paul Éluard