Description Usage Format Author(s) See Also Examples
Airborne counts of Fulmaris glacialis during the Aug/Sept 1998 and 1999 flights on the Dutch (Netherlands) part of the North Sea (NCP, Nederlands Continentaal Plat).
1 |
This data frame contains the following columns:
year of measurement: 1998 or 1999
x-coordinate in UTM zone 31
y-coordinate in UTM zone 31
sea water depth, in m
distance to coast of the Netherlands, in km.
observed density (number of birds per square km)
Dutch National Institute for Coastal and Marine Management (RIKZ), http://www.rikz.nl/
ncp.grid
E.J. Pebesma, R.N.M. Duin, P.A. Burrough, 2005. Mapping Sea Bird Densities over the North Sea: Spatially Aggregated Estimates and Temporal Changes. Environmetrics 16, (6), p 573-587.
1 2 3 4 5 6 |
year x y depth
Min. :1998 Min. :476210 Min. :5694947 Min. : 1.0
1st Qu.:1998 1st Qu.:535522 1st Qu.:5806777 1st Qu.:13.0
Median :1999 Median :568618 Median :5896021 Median :25.0
Mean :1999 Mean :576982 Mean :5889112 Mean :23.6
3rd Qu.:1999 3rd Qu.:604579 3rd Qu.:5946550 3rd Qu.:32.0
Max. :1999 Max. :739042 Max. :6150942 Max. :54.0
coast fulmar
Min. : 1.024 Min. : 0.000
1st Qu.: 4.857 1st Qu.: 0.000
Median : 38.696 Median : 0.000
Mean : 56.642 Mean : 1.005
3rd Qu.: 89.929 3rd Qu.: 0.000
Max. :263.205 Max. :46.487
demo(fulmar)
---- ~~~~~~
> # $Id: fulmar.R,v 1.3 2006-02-10 19:05:02 edzer Exp $
> library(sp)
> library(gstat)
> data(fulmar)
> data(ncp.grid)
> glm98 <- glm(formula = fulmar ~ depth + coast, family = quasipoisson,
+ data = fulmar[fulmar$year == 1998, ])
> glm99 <- glm(formula = fulmar ~ depth + coast, family = quasipoisson,
+ data = fulmar[fulmar$year == 1999, ])
> fulmar98 = data.frame(fulmar[fulmar$year == 1998,],
+ pr98 = predict(glm98, type = "response"))
> fulmar99 <- data.frame(fulmar[fulmar$year == 1999,],
+ pr99 = predict(glm99, type = "response"))
> pr98.grd <- predict(glm98, newdata = ncp.grid, type = "response", se.fit=TRUE)
> pr99.grd <- predict(glm99, newdata = ncp.grid, type = "response", se.fit=TRUE)
> pr <- data.frame(ncp.grid, pr98=pr98.grd$fit, pr99=pr99.grd$fit,
+ se98 = pr98.grd$se.fit, se99 = pr99.grd$se.fit)
> # B.3 create gstat object
> g <- gstat(id = "fulmar98", formula = fulmar~pr98, locations = ~x+y,
+ data = fulmar98, model = vgm(1.89629, "Exp", 50000, 0.852478),
+ beta = c(0,1), variance = "mu")
> g <- gstat(g, id = "fulmar99", formula = fulmar~pr99, locations = ~x+y,
+ data = fulmar99, model = vgm(2.52259, "Exp", 50000, 1.76474),
+ beta = c(0,1), variance = "mu")
> h <- g
> h <- gstat(h, id = c("fulmar98","fulmar99"),
+ model = vgm(2.18, "Exp", 50000, 1.22))
> # predict block means for blocks in ncp.grid$area (table 2; cokriging)
> library(maptools)
Checking rgeos availability: FALSE
Note: when rgeos is not available, polygon geometry computations in maptools depend on gpclib,
which has a restricted licence. It is disabled by default;
to enable gpclib, type gpclibPermit()
> areas.r = readShapePoly(system.file("external/ncp.shp", package="gstat"))
> #areas.r <- as.SpatialRings.Shapes(areas.shp$Shapes, areas.shp$att.data$WSVGEB_)
> coordinates(pr) = ~x+y
> #pr.df = overlay(pr, areas.r, fn = mean)
> pr.df = na.omit(as(aggregate(pr, areas.r, FUN = mean), "data.frame"))
> # match non-empty (and relevant) areas:
> #areas = SpatialPolygonsDataFrame(areas.r[c(2,3,4,16),"WSVGEB_"], pr.df[c(1,2,3,5),])#,match.ID=F)
> areas = SpatialPolygonsDataFrame(areas.r[c(1,2,12,7),"WSVGEB_"],
+ pr.df[c(1,2,3,5),], match.ID=F)
> # areas ID's 0 1 2 14
> sk = predict(g, areas)
[using simple kriging]
> cok = predict(h, areas)
Linear Model of Coregionalization found. Good.
[using simple cokriging]
> spplot(cok, c(3,5), names.attr = c("1998", "1999"),
+ main = "Fulmaris glacialis, density estimates\n(by irregular block cokriging)")
> sk = as.data.frame(sk)
> cok = as.data.frame(cok)
> print(data.frame(area = c(1,2,3,16),
+ SK98 = sk$fulmar98.pred, SE98 = sqrt(sk$fulmar98.var),
+ SK99 = sk$fulmar99.pred, SE99 = sqrt(sk$fulmar99.var),
+ CK98 = cok$fulmar98.pred, SE98 = sqrt(cok$fulmar98.var),
+ CK99 = cok$fulmar99.pred, SE99 = sqrt(cok$fulmar99.var)),
+ digits=3)
area SK98 SE98 SK99 SE99 CK98 SE98.1 CK99 SE99.1
1 1 2.9469 0.642 3.7043 0.740 2.9469 0.642 3.7043 0.740
2 2 0.2289 0.634 0.5242 0.731 0.2289 0.634 0.5242 0.731
3 3 0.0442 1.107 0.1003 1.276 0.0442 1.107 0.1003 1.276
4 16 0.0324 1.365 0.0517 1.574 0.0324 1.365 0.0517 1.574
> print(data.frame(area = c(1,2,3,16),
+ dSK = sk$fulmar99.pred - sk$fulmar98.pred,
+ SEdSK = sqrt(sk$fulmar98.var+sk$fulmar99.var),
+ dCOK = cok$fulmar99.pred - cok$fulmar98.pred,
+ SEdCOK = sqrt(cok$fulmar98.var+cok$fulmar99.var
+ - 2*cok$cov.fulmar98.fulmar99)),
+ digits=3)
area dSK SEdSK dCOK SEdCOK
1 1 0.7574 0.979 0.7574 0.113
2 2 0.2953 0.967 0.2953 0.112
3 3 0.0561 1.689 0.0561 0.195
4 16 0.0194 2.083 0.0194 0.240
There were 13 warnings (use warnings() to see them)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.