[R-sig-Geo] Density

Sat Mar 3 17:27:00 CET 2012

On Thu, 1 Mar 2012, alice schaffhauser wrote:

> Hi everybody, I have a data base on a plant species with records 
> attached to the centers of grid cells. Because I don't have softwares as 
> arcgis, I'd like to learn spatial analyses with R (and qgis?). I have 
> wgs84 coordinates and 3 status for the species population (casual=50 
> plants; unknown=500; established=1500).I just tried to see data with a 
> simple plot but I must represent the density of population.

It is impossible to understand your description of your input data. From 
your earlier post:

https://stat.ethz.ch/pipermail/r-sig-geo/2012-February/014291.html

to which you attached data (never a good idea, links to websites are 
always better), you could have done:

tab <- read.table("baseat_part.txt", header=TRUE, sep="\t")
library(sp)
coordinates(tab) <- c("WGS84_E", "WGS84_N")
proj4string(tab) <- CRS("+proj=longlat +datum=WGS84")
plot(tab)
spplot(tab, "status", col.regions=1:3)
table(tab$status)

However:

zerodist(tab)

reports many repeated points, so displaying the points by category is 
impossible to do sensibly, because they overlap.

> I found kde2d but I don't know if kernel suits and what it means. I 
> think that I can use it with the three populations status but I don't 
> know how.

No, you should not try methods unless you actually understand what they 
do. The function you refer to is a two-dimensional density plot, and is 
suited to point patterns, where one knows that all the objects belonging 
to the patterns have been observed. This looks like Austria, and assuming 
that it is, you could use tools for marked point patterns in spatstat, but 
preferably after projection, downloading national boundaries from GADM:

http://www.filefactory.com/file/c2a3441/n/AUT_adm0_RData

library(rgdal)
tab_lcc <- spTransform(tab, CRS("+init=epsg:3416"))
load("AUT_adm0.RData")
gadm_lcc <- spTransform(gadm, CRS("+init=epsg:3416"))
library(maptools)
library(spatstat)
spatstat.options(checkpolygons=FALSE)
# Polygon contains duplicated vertices and is self-intersecting
poly.owin <- as(gadm_lcc, "owin")
spatstat.options(checkpolygons=TRUE)
xy <- coordinates(tab_lcc)
points.ppp <- ppp(xy[,1], xy[,2], marks=tab_lcc$status, window=poly.owin)
# errors, some points not in Austria
plot(split(points.ppp), nrows=3)
Mask2 <- spsample(gadm_lcc, n=5000, type="regular")
XY <- coordinates(Mask2)
XY <- list(x=XY[,1], y=XY[,2])
Z <- density(split(points.ppp), xy=XY)
plot(Z, nrows=3)

Note that this is just an example, and only respects the point support of 
the data, not a grid size. The densities presuppose that all of the area 
has been surveyed with equal intensity in principle - this is a point 
pattern. If observations have not been made in any areas, they should be 
excluded from the window (here all of Austria).

Perhaps reading up on the methods involved would help?

Hope this clarifies,

Roger

> The grid resolution can vary ('1' means 3km of resolution and '2' means 
> 6 km so I must assign them to the first ones). I just heard about 
> krigeage (autokrige?) and nearest neighbors but I don't know how to use 
> them to calculate the number of plants by grid cell. The objective is to 
> have 'net cbf' format because I must transfer the distribution map to a 
> pollen emission map for a transport model. Thanks for your help! Alice 
> ------------------------------------------------ Alice SCHAFFHAUSER
>
> ------------------------------------------------
>
> tab <- read.table("base_at_part.txt", header=T, sep="\t")
> tab
> dim(tab)
> names(tab)
> summary(tab)
> plot(tab$WGS84_E,tab$WGS84_N, col="blue", xlab="WGS84 E", ylab="WGS84 N")
> ?
> #density histogram years
> hist(tab$year_record,?
> ???? probability=TRUE, breaks=20, col="light blue")
> points(density(tab$year_record, bw=1),? type='l', lwd=3, col='black')
> ?
> #density (kernel density estimation)
> library(MASS)?
> kde2d(tab$WGS84_E,tab$WGS84_N) -> tmp?
> filled.contour(tmp, color.palette = topo.colors)
> tmp
>
> ------------------------------------------------
> grid_cell??? status??? resolution??? year_record??? Y-LEFT??? X-LEFT??? y1??? y2??? x1??? x2??? WGS84_N??? WGS84_E??? x??? y
> 7764/2??? unknown??? 1??? 1999??? 4815??? 1625??? 48??? 15.5??? 16??? 25.5??? 48.2583333333333??? 16.425??? 16.4583333333??? 48.275
> 7658/1??? unknown??? 1??? 2005??? 4821??? 1520??? 48??? 21.5??? 15??? 20.5??? 48.3583333333333??? 15.3416666666667??? 15.375??? 48.375
> 7646/3??? unknown??? 1??? 2001??? 4818??? 1320??? 48??? 18.5??? 13??? 20.5??? 48.3083333333333??? 13.3416666666667??? 13.375??? 48.325
> 7850/2??? casual??? 1??? 2000??? 4809??? 1405??? 48??? 9.5??? 14??? 5.5??? 48.1583333333333??? 14.0916666666667??? 14.125??? 48.175
> 7745/1??? casual??? 1??? 1999??? 4815??? 1310??? 48??? 15.5??? 13??? 10.5??? 48.2583333333333??? 13.175??? 13.2083333333??? 48.275
> 7646/1??? unknown??? 1??? 2001??? 4821??? 1320??? 48??? 21.5??? 13??? 20.5??? 48.3583333333333??? 13.3416666666667??? 13.375??? 48.375
> 7842/4??? unknown??? 1??? 2002??? 4806??? 1245??? 48??? 6.5??? 12??? 45.5??? 48.1083333333333??? 12.7583333333333??? 12.7916666667??? 48.125
> 7446/2??? casual??? 2 ?? 2002??? 4833??? 1325??? 48??? 33.5??? 13??? 25.5??? 48.5583333333333??? 13.425??? 13.4583333333??? 48.575
> 7846/1??? casual??? 1??? 2002??? 4809??? 1320??? 48??? 9.5??? 13??? 20.5??? 48.1583333333333??? 13.3416666666667??? 13.375??? 48.175
> 7849/2??? unknown??? 1??? 2001??? 4809??? 1355??? 48??? 9.5??? 13??? 55.5??? 48.1583333333333??? 13.925??? 13.9583333333??? 48.175
> 7943/1??? unknown??? 1??? 2002??? 4803??? 1250??? 48??? 3.5??? 12??? 50.5??? 48.0583333333333??? 12.8416666666667??? 12.875??? 48.075
> 7646/4??? casual??? 1??? 2002??? 4818??? 1325??? 48??? 18.5??? 13??? 25.5??? 48.3083333333333??? 13.425??? 13.4583333333??? 48.325
> 8764/1??? unknown??? 2 ? 1974??? 4715??? 1620??? 47??? 15.5??? 16??? 20.5??? 47.2583333333333??? 16.3416666666667??? 16.375??? 47.275
> 7746/4??? established??? 1??? 2000??? 4812??? 1325??? 48??? 12.5??? 13??? 25.5??? 48.2083333333333??? 13.425??? 13.4583333333??? 48.225
> 7646/3??? unknown??? 1??? 2002??? 4818??? 1320??? 48??? 18.5??? 13??? 20.5??? 48.3083333333333??? 13.3416666666667??? 13.375??? 48.325
> 7745/4??? established??? 1??? 2002??? 4812??? 1315??? 48??? 12.5??? 13??? 15.5??? 48.2083333333333??? 13.2583333333333??? 13.2916666667??? 48.225
> 7747/2??? casual??? 1??? 2003??? 4815??? 1335??? 48??? 15.5??? 13??? 35.5??? 48.2583333333333??? 13.5916666666667??? 13.625??? 48.275
> 8146/2??? unknown??? 1??? 2003??? 4751??? 1325??? 47??? 51.5??? 13??? 25.5??? 47.8583333333333??? 13.425??? 13.4583333333??? 47.875
> 7646/1??? unknown??? 1??? 2003??? 4821??? 1320??? 48??? 21.5??? 13??? 20.5??? 48.3583333333333??? 13.3416666666667??? 13.375??? 48.375
> 7842/2??? unknown??? 1??? 2003??? 4809??? 1245??? 48??? 9.5??? 12??? 45.5??? 48.1583333333333??? 12.7583333333333??? 12.7916666667??? 48.175
> 7546/3??? unknown??? 1??? 2003??? 4824??? 1320??? 48??? 24.5??? 13??? 20.5??? 48.4083333333333??? 13.3416666666667??? 13.375??? 48.425
> 7752/4??? established??? 1??? 2003??? 4812??? 1425??? 48??? 12.5??? 14??? 25.5??? 48.2083333333333??? 14.425??? 14.4583333333??? 48.225
> 7863/4??? unknown??? 1??? 2004??? 4806??? 1615??? 48??? 6.5??? 16??? 15.5??? 48.1083333333333??? 16.2583333333333??? 16.2916666667??? 48.125
> 7963/1??? established??? 2 ?? 2003??? 4803??? 1610??? 48??? 3.5??? 16??? 10.5??? 48.0583333333333??? 16.175??? 16.2083333333??? 48.075
> 7849/2??? unknown??? 1??? 2003??? 4809??? 1355??? 48??? 9.5??? 13??? 55.5??? 48.1583333333333??? 13.925??? 13.9583333333??? 48.175
> 7752/4??? established??? 1??? 2003??? 4812??? 1425??? 48??? 12.5??? 14??? 25.5??? 48.2083333333333??? 14.425??? 14.4583333333??? 48.225
> 7850/4??? unknown??? 1??? 2003??? 4806??? 1405??? 48??? 6.5??? 14??? 5.5??? 48.1083333333333??? 14.0916666666667??? 14.125??? 48.125
> 7848/1??? unknown??? 2 ?? 2003??? 4809??? 1340??? 48??? 9.5??? 13??? 40.5??? 48.1583333333333??? 13.675??? 13.7083333333??? 48.175
> 7849/2??? unknown??? 1??? 2003??? 4809??? 1355??? 48??? 9.5??? 13??? 55.5??? 48.1583333333333??? 13.925??? 13.9583333333??? 48.175
> 	[[alternative HTML version deleted]]
>
>

-- 
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no