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

Andrew Vitale vitale232 at gmail.com
Fri Sep 12 01:01:35 CEST 2014

```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

-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)

--
*Andrew P. Vitale*
Masters Student
Department of Geography