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

[Andrew Vitale]

Thanks for the answers!  Robert's answer gives a tremendous increase in
speed.

Original code:
   user  system elapsed
  5.344   0.123   5.488

Robert's code:
   user  system elapsed
  0.170   0.002   0.173



On Thu, Sep 11, 2014 at 5:11 PM, Mathieu Basille <
basille.web at ase-research.org> wrote:

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


-- 
*Andrew P. Vitale*
Masters Student
Department of Geography
University of Nevada, Reno
vitale232 at gmail.com

	[[alternative HTML version deleted]]