[R-sig-Geo] Extracting informations from overlay

[Edzer Pebesma]

On 12/06/2011 05:17 PM, sandra werb wrote:
> 
> 
> 2011/12/6 Edzer Pebesma <edzer.pebesma at uni-muenster.de
> <mailto:edzer.pebesma at uni-muenster.de>>
> 
> 
> 
>     On 12/06/2011 12:06 PM, sandra werb wrote:
>     > Dear all,
>     >
>     > I am working with a Lidar dataset, which contains thousands of
>     > undistributed points. What I want to do is to count the number of
>     points
>     > within a raster cell and extract this information. Not just as single
>     > numbers, but together with the ID or the coordinates of the
>     raster. I came
>     > to this so far:
>     >
>     > - Create an empty raster, all cells are 0
>     >
>     > empty.r <- raster(nrows=733, ncols=732, xmn=-158, xmx=-121, ymn=-28.3,
>     > ymx=8.35, crs=NA)
>     > empty.r[] <- 0
>     >
>     > - read in the dataset of points containing only xy information and
>     define
>     > these as spatialpoints
>     >
>     > r1_points <- read.table("09_16_B.txt", sep="\t", dec=".", header=F)
>     > p <- r1_points[,-3]
>     >  p<-SpatialPoints(data.frame(p))
>     >
>     > - create an overlay of the raster and the spatial points and use
>     function
>     > table to
>     >
>     > o <- overlay(as(empty.r, "SpatialGrid"),p)
>     > tab <- table(o)
> 
>     untried, from ?over:
> 
>     g = as(empty.r, "SpatialGrid")
>     o = sapply(over(as(g, "SpatialPolygons"), p, returnList = TRUE), length)
>     gdf = SpatialGridDataFrame(g, data.frame(o = o))
>     spplot(gdf)
> 
> 
> Thanks, this works. When I sum up the counts within the overlay I get a
> number that is higher than the points I have in the dataset. Might it be
> the case that points, which lie on/very close to boundaries of the cells
> are counted twice?

Yes, that must be the case. (Maybe shifting all your points with 1 mm or
so helps?)

>  
> 
> 
>     >
>     > tab gives me the numbers of points within one cell, but I don't
>     know how to
>     > get this information together with the ID or XY coordinates as one
>     table or
>     > else. Can anyone help?
>     >
>     > Thank you very much in advance!
>     >
>     > Leni
>     >
>     >       [[alternative HTML version deleted]]
>     >
>     > _______________________________________________
>     > R-sig-Geo mailing list
>     > R-sig-Geo at r-project.org <mailto:R-sig-Geo at r-project.org>
>     > 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 <tel:%2B49%20251%208339763>
>      http://ifgi.uni-muenster.de
>     http://www.52north.org/geostatistics      e.pebesma at wwu.de
>     <mailto:e.pebesma at wwu.de>
> 
>     _______________________________________________
>     R-sig-Geo mailing list
>     R-sig-Geo at r-project.org <mailto:R-sig-Geo at r-project.org>
>     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      e.pebesma at wwu.de