Efficient way to obtain gridded count of overlapping polygons?

classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

Efficient way to obtain gridded count of overlapping polygons?

Lyndon Estes
Hello,

I am interested in creating a raster that summarizes the number of
overlapping polygons that underlie each cell (in order to display how
many species occur at that point).

I have come up with a function that seems to work, but I was wondering
if I had missed one that already exists, or if someone can recommend a
better way?

Herewith the code:


CountIntersections <- function(x, y, CID) {
# Counts number of polygons intersecting spatial points, specifically
designed for situations in which
# SpatialPolygons object contains overlapping polygons
# Args:
#   x: Input SpatialPointsDataFrame
#   y: SpatialPolygons* on which to count intersects
#   CID: Name of column containing unique polygon IDs
# Returns:
#   SpatialPointsDataFrame with column providing summary of number of
polygons with unique IDs intersected

  # Vector of unique IDs
  ids <- as.character(unique(y@data[, CID]))

  # lapply function in which each unique polygons is extracted and
overlaid by points, with 1 assigned for sp
  # presence, 0 for absence
  p.int <- lapply(ids, function(z) {
    p.ex <- y[y@data[, CID] == z, ]  # Extract polygon(s)
corresponding to unique ID extracted by lapply
    ov <- overlay(x, p.ex)  # Intersect points with sp polygon
    ov2 <- ov / ov  # Reduce intersected polygon IDs to 1
    ov2[is.na(ov2)] <- 0  # Convert NAs to 0
    return(ov2)
  })
  x$P.cnt <- Reduce("+", p.int)  # Sum list of intersect for each
species, and add value to dataframe of x
  return(x)
}

# Load spatial points (pts) and polygons of 3 species' ranges (sppX3)
extracted from IUCN mammals database
# downloaded here: http://www.iucnredlist.org/technical-documents/spatial-data
con <- url("http://sites.google.com/site/ldemisc/polygon-intersect/spp.shp.Rdata")
load(con)
close(con)

# Run function and convert to raster
pts.ct <- CountIntersections(pts, sppX3, "BINOMIAL")
pts.grd <- SpatialPixelsDataFrame(pts.ct, pts.ct@data)
pts.rst <- raster(pts.grd, layer = "P.cnt")

# Seems to work based on plots
plot(pts.rst)
plot(sppX3, add = T)

Thanks in advance for any advice.

Cheers, Lyndon

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Efficient way to obtain gridded count of overlapping polygons?

edzer
Yes, the (relatively) new method "over" does this for you, as in

pts$n = sapply(over(pts, geometry(sppX3), returnList = TRUE), length)

In the context of a running script:

require(sp)
con <- url(
"http://sites.google.com/site/ldemisc/polygon-intersect/spp.shp.Rdata")
load(con)
close(con)
proj4string(pts) = proj4string(sppX3) # otherwise over will fail here;
pts$n = sapply(over(pts, geometry(sppX3), returnList = TRUE), length)
gridded(pts)=TRUE
spplot(pts["n"], sp.layout=list("sp.polygons", sppX3, first=F))


On 02/17/2011 04:58 AM, Lyndon Estes wrote:

> Hello,
>
> I am interested in creating a raster that summarizes the number of
> overlapping polygons that underlie each cell (in order to display how
> many species occur at that point).
>
> I have come up with a function that seems to work, but I was wondering
> if I had missed one that already exists, or if someone can recommend a
> better way?
>
> Herewith the code:
>
>
> CountIntersections <- function(x, y, CID) {
> # Counts number of polygons intersecting spatial points, specifically
> designed for situations in which
> # SpatialPolygons object contains overlapping polygons
> # Args:
> #   x: Input SpatialPointsDataFrame
> #   y: SpatialPolygons* on which to count intersects
> #   CID: Name of column containing unique polygon IDs
> # Returns:
> #   SpatialPointsDataFrame with column providing summary of number of
> polygons with unique IDs intersected
>
>   # Vector of unique IDs
>   ids <- as.character(unique(y@data[, CID]))
>
>   # lapply function in which each unique polygons is extracted and
> overlaid by points, with 1 assigned for sp
>   # presence, 0 for absence
>   p.int <- lapply(ids, function(z) {
>     p.ex <- y[y@data[, CID] == z, ]  # Extract polygon(s)
> corresponding to unique ID extracted by lapply
>     ov <- overlay(x, p.ex)  # Intersect points with sp polygon
>     ov2 <- ov / ov  # Reduce intersected polygon IDs to 1
>     ov2[is.na(ov2)] <- 0  # Convert NAs to 0
>     return(ov2)
>   })
>   x$P.cnt <- Reduce("+", p.int)  # Sum list of intersect for each
> species, and add value to dataframe of x
>   return(x)
> }
>
> # Load spatial points (pts) and polygons of 3 species' ranges (sppX3)
> extracted from IUCN mammals database
> # downloaded here: http://www.iucnredlist.org/technical-documents/spatial-data
> con <- url("http://sites.google.com/site/ldemisc/polygon-intersect/spp.shp.Rdata")
> load(con)
> close(con)
>
> # Run function and convert to raster
> pts.ct <- CountIntersections(pts, sppX3, "BINOMIAL")
> pts.grd <- SpatialPixelsDataFrame(pts.ct, pts.ct@data)
> pts.rst <- raster(pts.grd, layer = "P.cnt")
>
> # Seems to work based on plots
> plot(pts.rst)
> plot(sppX3, add = T)
>
> Thanks in advance for any advice.
>
> Cheers, Lyndon
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> 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      [hidden email]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Efficient way to obtain gridded count of overlapping polygons?

Lyndon Estes
Hi Edzer,

Many thanks for the pointer.

My rough function is an order of magnitude slower than the 'over' one
you showed.

Sys.time() difference for my approach = 3.65 seconds, for yours =
0.323 seconds.

That will make a big difference when I run it for several thousand species.

Thanks again, Lyndon


On Thu, Feb 17, 2011 at 2:07 AM, Edzer Pebesma
<[hidden email]> wrote:

> Yes, the (relatively) new method "over" does this for you, as in
>
> pts$n = sapply(over(pts, geometry(sppX3), returnList = TRUE), length)
>
> In the context of a running script:
>
> require(sp)
> con <- url(
> "http://sites.google.com/site/ldemisc/polygon-intersect/spp.shp.Rdata")
> load(con)
> close(con)
> proj4string(pts) = proj4string(sppX3) # otherwise over will fail here;
> pts$n = sapply(over(pts, geometry(sppX3), returnList = TRUE), length)
> gridded(pts)=TRUE
> spplot(pts["n"], sp.layout=list("sp.polygons", sppX3, first=F))
>
>
> On 02/17/2011 04:58 AM, Lyndon Estes wrote:
>> Hello,
>>
>> I am interested in creating a raster that summarizes the number of
>> overlapping polygons that underlie each cell (in order to display how
>> many species occur at that point).
>>
>> I have come up with a function that seems to work, but I was wondering
>> if I had missed one that already exists, or if someone can recommend a
>> better way?
>>
>> Herewith the code:
>>
>>
>> CountIntersections <- function(x, y, CID) {
>> # Counts number of polygons intersecting spatial points, specifically
>> designed for situations in which
>> # SpatialPolygons object contains overlapping polygons
>> # Args:
>> #   x: Input SpatialPointsDataFrame
>> #   y: SpatialPolygons* on which to count intersects
>> #   CID: Name of column containing unique polygon IDs
>> # Returns:
>> #   SpatialPointsDataFrame with column providing summary of number of
>> polygons with unique IDs intersected
>>
>>   # Vector of unique IDs
>>   ids <- as.character(unique(y@data[, CID]))
>>
>>   # lapply function in which each unique polygons is extracted and
>> overlaid by points, with 1 assigned for sp
>>   # presence, 0 for absence
>>   p.int <- lapply(ids, function(z) {
>>     p.ex <- y[y@data[, CID] == z, ]  # Extract polygon(s)
>> corresponding to unique ID extracted by lapply
>>     ov <- overlay(x, p.ex)  # Intersect points with sp polygon
>>     ov2 <- ov / ov  # Reduce intersected polygon IDs to 1
>>     ov2[is.na(ov2)] <- 0  # Convert NAs to 0
>>     return(ov2)
>>   })
>>   x$P.cnt <- Reduce("+", p.int)  # Sum list of intersect for each
>> species, and add value to dataframe of x
>>   return(x)
>> }
>>
>> # Load spatial points (pts) and polygons of 3 species' ranges (sppX3)
>> extracted from IUCN mammals database
>> # downloaded here: http://www.iucnredlist.org/technical-documents/spatial-data
>> con <- url("http://sites.google.com/site/ldemisc/polygon-intersect/spp.shp.Rdata")
>> load(con)
>> close(con)
>>
>> # Run function and convert to raster
>> pts.ct <- CountIntersections(pts, sppX3, "BINOMIAL")
>> pts.grd <- SpatialPixelsDataFrame(pts.ct, pts.ct@data)
>> pts.rst <- raster(pts.grd, layer = "P.cnt")
>>
>> # Seems to work based on plots
>> plot(pts.rst)
>> plot(sppX3, add = T)
>>
>> Thanks in advance for any advice.
>>
>> Cheers, Lyndon
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> 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      [hidden email]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>



--
Lyndon Estes
Research Associate
Woodrow Wilson School
Princeton University
+1-609-258-2392 (o)
+1-609-258-6082 (f)
+1-202-431-0496 (m)
[hidden email]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Efficient way to obtain gridded count of overlapping polygons?

Robert Hijmans
In reply to this post by Lyndon Estes
> I am interested in creating a raster that summarizes the number of
> overlapping polygons that underlie each cell (in order to display how
> many species occur at that point).

Lyndon,

You can also do the below, with raster, see ?rasterize:

r <- rasterize(pols, r, fun='count')

or

r <- rasterize(pols, r, fun=function(x,...)length(na.omit(x)))

Reply | Threaded
Open this post in threaded view
|

Re: Efficient way to obtain gridded count of overlapping polygons?

Lyndon Estes
Hi Robert,

Many thanks for that suggestion as well, I appreciate it, which also
does the trick.

Just to update the time differences between the 3 methods, in case
anyone else is interested:

Function I wrote: 2.07 seconds
The "over" method from Edzer: 0.612 seconds
Rasterize function: 5.61 seconds (line of code used is below)


r2 <- rasterize(sppX3, r, fun = "count")

r was created from a shape of Africa with 1 degree resolution

Thanks again, Lyndon

On Thu, Feb 17, 2011 at 11:48 AM, Robert Hijmans <[hidden email]> wrote:

>
>> I am interested in creating a raster that summarizes the number of
>> overlapping polygons that underlie each cell (in order to display how
>> many species occur at that point).
>
> Lyndon,
>
> You can also do the below, with raster, see ?rasterize:
>
> r <- rasterize(pols, r, fun='count')
>
> or
>
> r <- rasterize(pols, r, fun=function(x,...)length(na.omit(x)))
>
>
> --
> View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Efficient-way-to-obtain-gridded-count-of-overlapping-polygons-tp6034590p6036670.html
> Sent from the R-sig-geo mailing list archive at Nabble.com.
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>



--
Lyndon Estes
Research Associate
Woodrow Wilson School
Princeton University
+1-609-258-2392 (o)
+1-609-258-6082 (f)
+1-202-431-0496 (m)
[hidden email]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo