All R code in this book

set.seed(1014)
options(digits = 3)

knitr::opts_chunk$set(
  comment = "#",
  collapse = knitr::is_latex_output(),
  cache = TRUE,
#  out.width = "70%",
  fig.align = 'center'
#  fig.width = 6,
#  fig.asp = 0.618,  # 1 / phi
#  fig.show = "hold"
)

options(dplyr.print_min = 6, dplyr.print_max = 6)
options(stars.crs = 17)
mapview::mapviewOptions(fgb = FALSE)

# for units: (?)
Sys.setenv(UDUNITS2_XML_PATH="")
if (knitr::is_latex_output())
    options(width = 66)
    #options(width = 72)
knitr::write_bib(c(
  "abind",
  "classInt",
  "colorspace",
  "gdalcubes",
  "gstat",
  "ggplot2",
  "ggspatial",
  "lwgeom",
  "osmar", 
  "raster",
  "rgee",
  "ecmwfr",
  "rnaturalearth",
  "RColorBrewer",
  "raster",
  "rstac",
  "sf", 
  "sp",
  "spacetime",
  "spatstat", #-->> generates too many authors!! FIXME: update before submit!
  "spdep",
  "splm",
  "stars",
  "stcos",
  "stplanr", 
  "terra",
  "tidyverse",
  "tsibble",
  "viridis",
  "units",
  "xts",
  "s2",
  "tmap",
  "mapview"
  ), "packages.bib", width = 60)
The first part of this book introduces concepts of spatial data science,
and uses R only to generate text output or figures. The R code used for
this is not shown, as it would distract from the message. The online
version of this book contains the R sections, which can be unfolded
on demand and copied into the clipboard for execution and experimenting.
The second part of this book explains how the concepts introduced in
part I are dealt with using R, and deals with basic handling and plotting
of spatial and spatiotemporal data. Part III is dedicated to statistical
modelling of spatial data.
This work is licensed under the
[Attribution-NonCommercial-NoDerivatives
4.0](https://creativecommons.org/licenses/by-nc-nd/4.0/legalcode)
International License.
set.seed(1014)
options(digits = 3)

knitr::opts_chunk$set(
  comment = "#",
  collapse = knitr::is_latex_output(),
  cache = TRUE,
#  out.width = "70%",
  fig.align = 'center'
#  fig.width = 6,
#  fig.asp = 0.618,  # 1 / phi
#  fig.show = "hold"
)

options(dplyr.print_min = 6, dplyr.print_max = 6)
options(stars.crs = 17)
mapview::mapviewOptions(fgb = FALSE)

# for units: (?)
Sys.setenv(UDUNITS2_XML_PATH="")
if (knitr::is_latex_output())
    options(width = 66)
    #options(width = 72)
Text introducing Part I
library(tidyverse)
library(sf)
system.file("gpkg/nc.gpkg", package="sf") %>%
    read_sf() -> nc
nc.32119 <- st_transform(nc, 32119)
nc.32119 %>%
    select(BIR74) %>%
    plot(graticule = TRUE, axes = TRUE)
nc %>% select(AREA, BIR74, SID74) %>% print(n = 3)
year_labels = c("SID74" = "1974 - 1978", "SID79" = "1979 - 1984")
nc.32119 %>% select(SID74, SID79) %>% gather(VAR, SID, -geom) -> nc2
ggplot() + geom_sf(data = nc2, aes(fill = SID)) + 
  facet_wrap(~VAR, ncol = 1, labeller = labeller(VAR = year_labels)) +
  scale_y_continuous(breaks = 34:36) +
  scale_fill_gradientn(colors = sf.colors(20)) +
  theme(panel.grid.major = element_line(color = "white"))
suppressPackageStartupMessages(library(mapview))
nc.32119 %>% mapview(zcol = "BIR74", legend = TRUE, col.regions = sf.colors)
library(stars)
par(mfrow = c(2, 2))
par(mar = rep(1, 4))
tif <- system.file("tif/L7_ETMs.tif", package = "stars")
x <- read_stars(tif)[,,,1]
image(x, main = "(a)")
image(x[,1:10,1:10], text_values = TRUE, border = 'grey', main = "(b)")
image(x, main = "(c)")
set.seed(131)
pts = st_sample(st_as_sfc(st_bbox(x)), 3)
plot(pts, add = TRUE, pch = 3, col = 'blue')
image(x, main = "(d)")
plot(st_buffer(pts, 500), add = TRUE, pch = 3, border = 'blue', col = NA, lwd = 2)
st_extract(x, pts)
aggregate(x, st_buffer(pts, 500), FUN = mean) %>% st_as_sf()
plot(st_rasterize(nc["BIR74"], dx = 0.1), col =sf.colors(), breaks = "equal")
x = 1:5
y = 1:4
d = st_dimensions(x = x, y = y, .raster = c("x", "y"))
m = matrix(runif(20),5,4)
r1 = st_as_stars(r = m, dimensions = d)

r = attr(d, "raster")
r$affine = c(0.2, -0.2)
attr(d, "raster") = r
r2 = st_as_stars(r = m, dimensions = d)

r = attr(d, "raster")
r$affine = c(0.1, -0.3)
attr(d, "raster") = r
r3 = st_as_stars(r = m, dimensions = d)

x = c(1, 2, 3.5, 5, 6)
y = c(1, 1.5, 3, 3.5)
d = st_dimensions(x = x, y = y, .raster = c("x", "y"))
r4 = st_as_stars(r = m, dimensions = d)

grd = st_make_grid(cellsize = c(10,10), offset = c(-130,10), n= c(8,5), crs=st_crs(4326))
r5 = st_transform(grd, "+proj=laea +lon_0=-70 +lat_0=35")

par(mfrow = c(2,3), mar = c(0.1, 1, 1.1, 1))
r1 = st_make_grid(cellsize = c(1,1), n = c(5,4), offset = c(0,0))
plot(r1, main = "regular")
plot(st_geometry(st_as_sf(r2)), main = "rotated")
plot(st_geometry(st_as_sf(r3)), main = "sheared")
plot(st_geometry(st_as_sf(r4, as_points = FALSE)), main = "rectilinear")
plot(st_geometry((r5)), main = "curvilinear")
knitr::include_graphics("images/sf_deps.png")
set.seed(1014)
options(digits = 3)

knitr::opts_chunk$set(
  comment = "#",
  collapse = knitr::is_latex_output(),
  cache = TRUE,
#  out.width = "70%",
  fig.align = 'center'
#  fig.width = 6,
#  fig.asp = 0.618,  # 1 / phi
#  fig.show = "hold"
)

options(dplyr.print_min = 6, dplyr.print_max = 6)
options(stars.crs = 17)
mapview::mapviewOptions(fgb = FALSE)

# for units: (?)
Sys.setenv(UDUNITS2_XML_PATH="")
if (knitr::is_latex_output())
    options(width = 66)
    #options(width = 72)
par(mar = rep(0,4))
plot(3, 4, xlim = c(-6,6), ylim = c(-6,6), asp = 1)
axis(1, pos = 0, at = 0:6)
axis(2, pos = 0, at = -6:6)
xd = seq(-5, 5, by = .1)
lines(xd, sqrt(25 - xd^2), col = 'grey')
lines(xd, -sqrt(25 - xd^2), col = 'grey')
arrows(0, 0, 3, 4, col = 'red', length = .15, angle = 20)
text(1.5, 2.7, label = "r", col = 'red')
xd = seq(3/5, 1, by = .1)
lines(xd, sqrt(1 - xd^2), col = 'red')
text(1.2, 0.5, label = parse(text = "phi"), col = 'red')
lines(c(3,3), c(0,4), lty = 2, col = 'blue')
lines(c(0,3), c(4,4), lty = 2, col = 'blue')
text(3.3, 0.3, label = "x", col = 'blue')
text(0.3, 4.3, label = "y", col = 'blue')
suppressPackageStartupMessages(library(sf))
e = cbind(-90:90,0) # equator
f1 = rbind(cbind(0, -90:90)) # 0/antimerid.
f2 = rbind(cbind(90, -90:90), cbind(270, 90:-90))# +/- 90
eq = st_sfc(st_linestring(e), st_linestring(f1), st_linestring(f2), crs=4326)

geoc = st_transform(eq, "+proj=geocent")
cc = rbind(geoc[[1]], NA, geoc[[2]], NA, geoc[[3]])
from3d = function(x, offset, maxz, minz) {
    x = x[,c(2,3,1)] + offset # move to y right, x up, z backw
    x[,2] = x[,2] - maxz      # shift y to left
    d = maxz
    z = x[,3] - minz + offset
    x[,1] = x[,1] * (d/z)
    x[,2] = x[,2] * (d/z)
    x[,1:2]
}
maxz = max(cc[,3], na.rm = TRUE)
minz = min(cc[,3], na.rm = TRUE)
offset = 3e7
circ = from3d(cc, offset, maxz, minz)
mx = max(cc, na.rm = TRUE) * 1.1
x = rbind(c(0, 0, 0), c(mx, 0, 0))
y = rbind(c(0, 0, 0), c(0, mx, 0))
z = rbind(c(0, 0, 0), c(0, 0, mx))
ll = rbind(x, NA, y, NA, z)
l0 =  from3d(ll, offset, maxz, minz)
mx = max(cc, na.rm = TRUE) * 1.2
x = rbind(c(0, 0, 0), c(mx, 0, 0))
y = rbind(c(0, 0, 0), c(0, mx, 0))
z = rbind(c(0, 0, 0), c(0, 0, mx))
ll = rbind(x, NA, y, NA, z)
l =  from3d(ll, offset, maxz, minz)

par(mfrow = c(1, 2))
par(mar=rep(0,4))
plot.new()
plot.window(xlim = c(min(circ[,1],na.rm = TRUE), 3607103*1.02), 
                        ylim = c(min(circ[,2],na.rm = TRUE), 2873898*1.1), asp = 1)
lines(circ)
lines(l0)
text(l[c(2,5,8),], labels = c("x", "y", "z"), col = 'red')
# add POINT(60 47)
p = st_as_sfc("POINT(60 47)", crs = 4326) %>% st_transform("+proj=geocent")
p = p[[1]]
pts = rbind(c(0,0,0), c(p[1],0,0), c(p[1],p[2],0), c(p[1],p[2],p[2]))
ptsl =  from3d(pts, offset, maxz, minz)
lines(ptsl, col = 'blue', lty = 2, lwd = 2)
points(ptsl[4,1], ptsl[4,2], col = 'blue', cex = 1, pch = 16)

plot.new()
plot.window(xlim = c(min(circ[,1],na.rm = TRUE), 3607103*1.02), 
                        ylim = c(min(circ[,2],na.rm = TRUE), 2873898*1.1), asp = 1)
lines(circ)

p = st_as_sfc("POINT(60 47)", crs = 4326) %>% st_transform("+proj=geocent")
p = p[[1]]
pts = rbind(c(0,0,0), c(p[1],p[2],p[3]))
pt =  from3d(pts, offset, maxz, minz)
lines(pt)
points(pt[2,1], pt[2,2], col = 'blue', cex = 1, pch = 16)

p0 = st_as_sfc("POINT(60 0)", crs = 4326) %>% st_transform("+proj=geocent")
p0 = p0[[1]]
pts = rbind(c(0,0,0), c(p0[1],p0[2],p0[3]))
pt =  from3d(pts, offset, maxz, minz)
lines(pt)

p0 = st_as_sfc("POINT(0 0)", crs = 4326) %>% st_transform("+proj=geocent")
p0 = p0[[1]]
pts = rbind(c(0,0,0), c(p0[1],p0[2],p0[3]))
pt =  from3d(pts, offset, maxz, minz)
lines(pt)

p0 = st_as_sfc("POINT(0 90)", crs = 4326) %>% st_transform("+proj=geocent")
p0 = p0[[1]]
pts = rbind(c(0,0,0), c(p0[1],p0[2],p0[3]))
pt =  from3d(pts, offset, maxz, minz)
lines(pt, lty = 2)

p0 = st_as_sfc("POINT(90 0)", crs = 4326) %>% st_transform("+proj=geocent")
p0 = p0[[1]]
pts = rbind(c(0,0,0), c(p0[1],p0[2],p0[3]))
pt =  from3d(pts, offset, maxz, minz)
lines(pt, lty = 2)

f1 = rbind(cbind(0:60, 0))
arc = st_sfc(st_linestring(f1), crs=4326)
geoc = st_transform(arc, "+proj=geocent")
cc = rbind(geoc[[1]])
circ = from3d(cc, offset, maxz, minz)
lines(circ, col = 'red', lwd = 2, lty = 2)

f1 = rbind(cbind(60, 0:47))
arc = st_sfc(st_linestring(f1), crs=4326)
geoc = st_transform(arc, "+proj=geocent")
cc = rbind(geoc[[1]])
circ = from3d(cc, offset, maxz, minz)
lines(circ, col = 'blue', lwd = 2, lty = 2)

text(pt[1,1]+100000, pt[1,2]+50000, labels = expression(phi), col = 'blue') # lat
text(pt[1,1]+20000, pt[1,2]-50000, labels = expression(lambda), col = 'red') # lng
p = st_as_sfc("POINT(60 47)", crs = 4326)
p[[1]]
p = st_as_sfc("POINT(60 47)", crs = 4326) %>% st_transform("+proj=geocent")
p[[1]]
par(mar = rep(0,4))
x = 4
y = 5/8 * sqrt(48)
plot(x, y, xlim = c(-6,6), ylim = c(-8,8), asp = 1)
axis(1, pos = 0, at = 0:9)
axis(2, pos = 0, at = -5:5)
xd = seq(-8, 8, by = .1)
lines(xd, 5/8 * sqrt(64 - xd^2), col = 'grey')
lines(xd, 5/8 * -sqrt(64 - xd^2), col = 'grey')
arrows(0, 0, x, y, col = 'red', length = .15, angle = 20)
b = (x * 25) / (-y * 64)
a = y - x * b
abline(a, b, col = 'grey')
b = -1/b
x0 = x - y / b
arrows(x0, 0, x, y, col = 'blue', length = .15, angle = 20)
text(1.2, 0.5, label = parse(text = "psi"), col = 'red')
text(3, 0.5, label = parse(text = "phi"), col = 'blue')
pts = st_sfc(st_point(c(13.4050, 52.5200)), st_point(c(2.3522, 48.8566)), crs = 'EPSG:4326')
sf_use_s2(FALSE)
d1 = c(gc_ellipse = st_distance(pts)[1,2])
sf_use_s2(TRUE)
# or, without using s2, use st_distance(st_transform(pts, "+proj=longlat +ellps=sphere"))
d2 = c(gc_sphere = st_distance(pts)[1,2])
p = st_transform(pts, "+proj=geocent")
d3 = c(straight_ellipse = units::set_units(sqrt(sum(apply(do.call(cbind, p), 1, diff)^2)), m))
p2 = st_transform(pts, "+proj=longlat +ellps=sphere") %>% st_transform("+proj=geocent")
d4 = c(straight_sphere = units::set_units(sqrt(sum(apply(do.call(cbind, p2), 1, diff)^2)), m))
res = c(d1,d3,d2,d4)
# print as km, re-add names:
res %>% units::set_units(km) %>% setNames(names(res)) %>% print(digits = 5)
library(stars)
library(rnaturalearth)
countries110 = st_as_sf(countries110)
uk = countries110[countries110$admin %in% c("United Kingdom"),] %>%
        st_geometry()
r = read_stars("data/uk_os_OSTN15_NTv2_OSGBtoETRS.tif")
# r = read_stars("/vsicurl/https://cdn.proj.org/uk_os_OSTN15_NTv2_OSGBtoETRS.tif")
hook = function() {
        plot(uk, border = "orange", col = NA, add = TRUE)
}
plot(r[,,,1:2], axes = TRUE, hook = hook, key.pos = 4)
h = read_stars("data/uk_os_OSGM15_GB.tif")
# h = read_stars("/vsicurl/https://cdn.proj.org/uk_os_OSGM15_GB.tif")
plot(h, axes = TRUE, reset = FALSE)
plot(uk, border = "orange", col = NA, add = TRUE)
set.seed(1014)
options(digits = 3)

knitr::opts_chunk$set(
  comment = "#",
  collapse = knitr::is_latex_output(),
  cache = TRUE,
#  out.width = "70%",
  fig.align = 'center'
#  fig.width = 6,
#  fig.asp = 0.618,  # 1 / phi
#  fig.show = "hold"
)

options(dplyr.print_min = 6, dplyr.print_max = 6)
options(stars.crs = 17)
mapview::mapviewOptions(fgb = FALSE)

# for units: (?)
Sys.setenv(UDUNITS2_XML_PATH="")
if (knitr::is_latex_output())
    options(width = 66)
    #options(width = 72)
library(sf)
par(mfrow = c(2,4))
par(mar = c(1,1,1.2,1))

# 1
p = st_point(0:1)
plot(p, pch = 16)
title("point")
box(col = 'grey')

# 2
mp = st_multipoint(rbind(c(1,1), c(2, 2), c(4, 1), c(2, 3), c(1,4)))
plot(mp, pch = 16)
title("multipoint")
box(col = 'grey')

# 3
ls = st_linestring(rbind(c(1,1), c(5,5), c(5, 6), c(4, 6), c(3, 4), c(2, 3)))
plot(ls, lwd = 2)
title("linestring")
box(col = 'grey')

# 4
mls = st_multilinestring(list(
  rbind(c(1,1), c(5,5), c(5, 6), c(4, 6), c(3, 4), c(2, 3)),
  rbind(c(3,0), c(4,1), c(2,1))))
plot(mls, lwd = 2)
title("multilinestring")
box(col = 'grey')

# 5 polygon
po = st_polygon(list(rbind(c(2,1), c(3,1), c(5,2), c(6,3), c(5,3), c(4,4), c(3,4), c(1,3), c(2,1)),
    rbind(c(2,2), c(3,3), c(4,3), c(4,2), c(2,2))))
plot(po, border = 'black', col = '#ff8888', lwd = 2)
title("polygon")
box(col = 'grey')

# 6 multipolygon
mpo = st_multipolygon(list(
    list(rbind(c(2,1), c(3,1), c(5,2), c(6,3), c(5,3), c(4,4), c(3,4), c(1,3), c(2,1)),
        rbind(c(2,2), c(3,3), c(4,3), c(4,2), c(2,2))),
    list(rbind(c(3,7), c(4,7), c(5,8), c(3,9), c(2,8), c(3,7)))))
plot(mpo, border = 'black', col = '#ff8888', lwd = 2)
title("multipolygon")
box(col = 'grey')

# 7 geometrycollection
gc = st_geometrycollection(list(po, ls + c(0,5), st_point(c(2,5)), st_point(c(5,4))))
plot(gc, border = 'black', col = '#ff6666', pch = 16, lwd = 2)
title("geometrycollection")
box(col = 'grey')
p
mp
ls
mls
po
mpo
gc
(ls = st_linestring(rbind(c(0,0), c(1,1), c(2,2), c(0,2), c(1,1), c(2,0))))
c(is_simple = st_is_simple(ls))
st_point(c(1,3,2))
st_point(c(1,3,2), dim = "XYM")
st_linestring(rbind(c(3,1,2,4), c(4,4,2,2)))
(e = st_intersection(st_point(c(0,0)), st_point(c(1,1))))
st_point()
st_linestring(matrix(1,1,3)[0,], dim = "XYM")
library(sf)
polygon = po = st_polygon(list(rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0))))
p0 = st_polygon(list(rbind(c(-1,-1), c(2,-1), c(2,2), c(-1,2), c(-1,-1))))
line = li = st_linestring(rbind(c(.5, -.5), c(.5, 0.5)))
s = st_sfc(po, li)

par(mfrow = c(3,3))
par(mar = c(1,1,1,1))

# "1020F1102"
# 1: 1
plot(s, col = c(NA, 'darkgreen'), border = 'blue', main = expression(paste("I(pol)",intersect(),"I(line) = 1")))
lines(rbind(c(.5,0), c(.5,.495)), col = 'red', lwd = 2)
points(0.5, 0.5, pch = 1)

# 2: 0
plot(s, col = c(NA, 'darkgreen'), border = 'blue', main = expression(paste("I(pol)",intersect(),"B(line) = 0")))
points(0.5, 0.5, col = 'red', pch = 16)

# 3: 2
plot(s, col = c(NA, 'darkgreen'), border = 'blue', main = expression(paste("I(pol)",intersect(),"E(line) = 2")))
plot(po, col = '#ff8888', add = TRUE)
plot(s, col = c(NA, 'darkgreen'), border = 'blue', add = TRUE)

# 4: 0
plot(s, col = c(NA, 'darkgreen'), border = 'blue', main = expression(paste("B(pol)",intersect(),"I(line) = 0")))
points(.5, 0, col = 'red', pch = 16)

# 5: F
plot(s, col = c(NA, 'darkgreen'), border = 'blue', main = expression(paste("B(pol)",intersect(),"B(line) = F")))

# 6: 1
plot(s, col = c(NA, 'darkgreen'), border = 'blue', main = expression(paste("B(pol)",intersect(),"E(line) = 1")))
plot(po, border = 'red', col = NA, add = TRUE, lwd = 2)

# 7: 1
plot(s, col = c(NA, 'darkgreen'), border = 'blue', main = expression(paste("E(pol)",intersect(),"I(line) = 1")))
lines(rbind(c(.5, -.5), c(.5, 0)), col = 'red', lwd = 2)

# 8: 0
plot(s, col = c(NA, 'darkgreen'), border = 'blue', main = expression(paste("E(pol)",intersect(),"B(line) = 0")))
points(.5, -.5, col = 'red', pch = 16)

# 9: 2
plot(s, col = c(NA, 'darkgreen'), border = 'blue', main = expression(paste("E(pol)",intersect(),"E(line) = 2")))
plot(p0 / po, col = '#ff8888', add = TRUE)
plot(s, col = c(NA, 'darkgreen'), border = 'blue', add = TRUE)
st_relate(polygon, line)
par(mar = rep(0,4), mfrow = c(1, 2))
plot(st_geometry(nc)[1], col = NA, border = 'black')
plot(st_convex_hull(st_geometry(nc)[1]), add = TRUE, col = NA, border = 'red')
box()
set.seed(131)
mp = st_multipoint(matrix(runif(20), 10))
plot(mp)
plot(st_voronoi(mp), add = TRUE, col = NA, border = 'red')
box()
par(mar = rep(.1, 4), mfrow = c(1, 2))
sq = function(pt, sz = 1) st_polygon(list(rbind(c(pt - sz), 
  c(pt[1] + sz, pt[2] - sz), c(pt + sz), c(pt[1] - sz, pt[2] + sz), c(pt - sz))))
x = st_sf(box = 1:3, st_sfc(sq(c(0,0)), sq(c(1.7, -0.5)), sq(c(0.5, 1))))
plot(st_geometry(x), col = NA, border = sf.colors(3, categorical=TRUE), lwd = 3)
plot(st_intersection(st_geometry(x)), col = sf.colors(7, categorical=TRUE, alpha = .5))
par(mar = rep(.1, 4), mfrow = c(1, 2)) 
xg = st_geometry(x)
plot(st_difference(xg), col = sf.colors(3, alpha = .5, categorical=TRUE))
plot(st_difference(xg[3:1]), col = sf.colors(3, alpha = .5, categorical=TRUE))
library(stars)
ls = st_sf(a = 1:2, st_sfc(st_linestring(rbind(c(0.1, 0), c(1, .9)))))
grd = st_as_stars(st_bbox(ls), nx = 10, ny = 10, xlim = c(0, 1.0), ylim = c(0, 1),
   values = -1)
r = st_rasterize(ls, grd, options = "ALL_TOUCHED=TRUE")
r[r == -1] = NA
plot(st_geometry(st_as_sf(grd)), border = 'orange', col = NA, 
     reset = FALSE, key.pos=NULL)
plot(r, axes = TRUE, add = TRUE) # ALL_TOUCHED=FALSE;
plot(ls, add = TRUE, col = "red", lwd = 2)
set.seed(1014)
options(digits = 3)

knitr::opts_chunk$set(
  comment = "#",
  collapse = knitr::is_latex_output(),
  cache = TRUE,
#  out.width = "70%",
  fig.align = 'center'
#  fig.width = 6,
#  fig.asp = 0.618,  # 1 / phi
#  fig.show = "hold"
)

options(dplyr.print_min = 6, dplyr.print_max = 6)
options(stars.crs = 17)
mapview::mapviewOptions(fgb = FALSE)

# for units: (?)
Sys.setenv(UDUNITS2_XML_PATH="")
if (knitr::is_latex_output())
    options(width = 66)
    #options(width = 72)
library(sf)
suppressPackageStartupMessages(library(maps))
# maps:
par(mfrow = c(2,2))
par(mar = c(1,1.2,1,1))
m = st_as_sf(map(fill=TRUE, plot=FALSE))
a = m[m$ID == "Antarctica", ]
st_bbox(a)
library(s2)
s2_bounds_cap(a)
s2_bounds_rect(a)
st_bbox(m[m$ID == "Fiji",])
s2_bounds_rect(m[m$ID == "Fiji",])
library(sf)
suppressPackageStartupMessages(library(maps))
# maps:
par(mfrow = c(2,2))
par(mar = c(1,1.2,1,1))
m = st_as_sf(map(fill=TRUE, plot=FALSE))
m = m[m$ID == "Antarctica", ]
plot(st_geometry(m), asp = 2)
title("a (not valid)")
# ne:
library(rnaturalearth)
ne = ne_countries(returnclass = "sf")
ne = ne[ne$region_un == "Antarctica", "region_un"]
plot(st_geometry(ne), asp = 2)
title("b (valid)")
# 3031
m %>%
  st_geometry() %>%
  st_transform(3031) %>%
  plot()
title("c (valid)")
ne %>%
  st_geometry() %>%
  st_transform(3031) %>%
  plot()
title("d (not valid)")
set.seed(1014)
options(digits = 3)

knitr::opts_chunk$set(
  comment = "#",
  collapse = knitr::is_latex_output(),
  cache = TRUE,
#  out.width = "70%",
  fig.align = 'center'
#  fig.width = 6,
#  fig.asp = 0.618,  # 1 / phi
#  fig.show = "hold"
)

options(dplyr.print_min = 6, dplyr.print_max = 6)
options(stars.crs = 17)
mapview::mapviewOptions(fgb = FALSE)

# for units: (?)
Sys.setenv(UDUNITS2_XML_PATH="")
if (knitr::is_latex_output())
    options(width = 66)
    #options(width = 72)
library(sf)
library(dplyr)
system.file("gpkg/nc.gpkg", package="sf") %>%
    read_sf() %>%
    st_transform(32119) %>%
    select(BIR74, SID74, NAME) %>%
    st_centroid() -> x
nc <- read_sf(system.file("gpkg/nc.gpkg", package="sf"))
#plot(nc["SID74"], axes = TRUE)
sf_use_s2(TRUE)
# encode quadrant by two logicals:
nc$lng = st_coordinates(st_centroid(st_geometry(nc)))[,1] > -79
nc$lat = st_coordinates(st_centroid(st_geometry(nc)))[,2] > 35.5
nc.grp = aggregate(nc["SID74"], list(nc$lng, nc$lat), sum)
plot(nc.grp["SID74"], axes = TRUE)
nc <- st_transform(nc, 2264)
gr = st_sf(
   label = apply(expand.grid(1:10, LETTERS[10:1])[,2:1], 1, paste0, collapse = " "),
   geom = st_make_grid(nc))
plot(st_geometry(nc), reset = FALSE, border = 'grey')
plot(st_geometry(gr), add = TRUE)
a = aggregate(nc["SID74"], gr, sum)
c(sid74_sum_counties = sum(nc$SID74), sid74_sum_rectangles = sum(a$SID74, na.rm = TRUE))
set.seed(1014)
options(digits = 3)

knitr::opts_chunk$set(
  comment = "#",
  collapse = knitr::is_latex_output(),
  cache = TRUE,
#  out.width = "70%",
  fig.align = 'center'
#  fig.width = 6,
#  fig.asp = 0.618,  # 1 / phi
#  fig.show = "hold"
)

options(dplyr.print_min = 6, dplyr.print_max = 6)
options(stars.crs = 17)
mapview::mapviewOptions(fgb = FALSE)

# for units: (?)
Sys.setenv(UDUNITS2_XML_PATH="")
if (knitr::is_latex_output())
    options(width = 66)
    #options(width = 72)
# (C) 2019, Edzer Pebesma, CC-BY-SA
set.seed(1331)
suppressPackageStartupMessages(library(stars))
library(colorspace)
tif = system.file("tif/L7_ETMs.tif", package = "stars")
r = read_stars(tif)

nrow = 5
ncol = 8
m = r[[1]][1:nrow,1:ncol,1]
dim(m) = c(x = nrow, y = ncol) # named dim
s = st_as_stars(m)
# s
attr(s, "dimensions")[[1]]$delta = 3
attr(s, "dimensions")[[2]]$delta = -.5
attr(attr(s, "dimensions"), "raster")$affine = c(-1.2, 0.0)

plt = function(x, yoffset = 0, add, li = TRUE) {
    attr(x, "dimensions")[[2]]$offset = attr(x, "dimensions")[[2]]$offset + yoffset
    l = st_as_sf(x, as_points = FALSE)
    pal = sf.colors(10)
    if (li)
        pal = lighten(pal, 0.3 + rnorm(1, 0, 0.1))
    if (! add)
        plot(l, axes = FALSE, breaks = "equal", pal = pal, reset = FALSE, border = grey(.75), key.pos = NULL, main = NULL, xlab = "time")
    else
        plot(l, axes = TRUE, breaks = "equal", pal = pal, add = TRUE, border = grey(.75))
    u = st_union(l)
    # print(u)
    plot(st_geometry(u), add = TRUE, col = NA, border = 'black', lwd = 2.5)
}

pl = function(s, x, y, add = TRUE, randomize = FALSE) {
  attr(s, "dimensions")[[1]]$offset = x
  attr(s, "dimensions")[[2]]$offset = y
  m = r[[1]][y + 1:nrow,x + 1:ncol,1]
  if (randomize)
    m = m[sample(y + 1:nrow),x + 1:ncol]
  dim(m) = c(x = nrow, y = ncol) # named dim
  s[[1]] = m
  plt(s, 0, add)
  plt(s, 1, TRUE)
  plt(s, 2, TRUE)
  plt(s, 3, TRUE)
  plt(s, 4, TRUE)
  plt(s, 5, TRUE)
  plt(s, 6, TRUE)
  plt(s, 7, TRUE)
  plt(s, 8, TRUE, FALSE)
}

plot.new()
par(mar = rep(0.5,4))
plot.window(xlim = c(-12,15), ylim = c(-5,10), asp=1)
pl(s, 0, 0)
# box()
text(-10, 0, "time", srt = -90, col = 'black')
text(-5,  6.5, "latitude", srt = 25, col = 'black')
text( 5,  8.5, "longitude", srt = 0, col = 'black')
# (C) 2021, Jonathan Bahlmann, CC-BY-SA
#  https://github.com/Open-EO/openeo.org/tree/datacube_doc/documentation/1.0/datacubes/.scripts
# based on work by Edzer Pebesma, 2019, here: https://gist.github.com/edzer/5f1b0faa3e93073784e01d5a4bb60eca

# plotting runs via a dummy stars object with x, y dimensions (no bands)
# to not be overly dependent on an input image, time steps and bands
# are displayed by replacing the matrix contained in the dummy stars object
# every time something is plotted

# packages, read input ----
set.seed(1331)
library(stars)
suppressPackageStartupMessages(library(colorspace))
suppressPackageStartupMessages(library(scales))

# make color palettes ----
blues <- sequential_hcl(n = 20, h1 = 211, c1 = 80, l1 = 40, l2 = 100, p1 = 2)
greens <- sequential_hcl(n = 20, h1 = 134, c1 = 80, l1 = 40, l2 = 100, p1 = 2)
reds <- sequential_hcl(n = 20, h1 = 360, c1 = 80, l1 = 40, l2 = 100, p1 = 2)
purples <- sequential_hcl(n = 20, h1 = 299, c1 = 80, l1 = 40, l2 = 100, p1 = 2)
greys <- sequential_hcl(n = 20, h1 = 0, c1 = 0, l1 = 40, l2 = 100, p1 = 2)

# matrices from raster ----
# make input matrices from an actual raster image
input <- read_stars("data/iceland_delta_cutout_2.tif") # this raster needs approx 6x7 format
# if the input raster is changed, every image where a pixel value is written as text needs to be checked and corrected accordingly
input <- input[,,,1:4]
warped <- st_warp(input, crs = st_crs(input), cellsize = 200) # warp to approx. 6x7 pixel

# these are only needed for resampling
warped_highres <- st_warp(warped, crs = st_crs(warped), cellsize = 100) # with different input, cellsize must be adapted
# this is a bit of a trick, because 3:4 is different format than 6:7
# when downsampling, the raster of origin isn't so important anyway
warped_lowres <- st_warp(warped_highres[,1:11,,], crs = st_crs(warped), cellsize = 390)
# plot(warped_lowres)
# image(warped[,,,1], text_values = TRUE)

t1 <- floor(matrix(runif(42, -30, 150), ncol = 7)) # create timesteps 2 and 3 randomly
t2 <- floor(matrix(runif(42, -250, 50), ncol = 7))

# create dummy stars object ----
make_dummy_stars <- function(x, y, d1, d2, aff) {
  m = warped_highres[[1]][1:x,1:y,1] # underlying raster doesn't matter because it's just dummy construct
  dim(m) = c(x = x, y = y) # named dim
  dummy = st_as_stars(m)
  attr(dummy, "dimensions")[[1]]$delta = d1
  attr(dummy, "dimensions")[[2]]$delta = d2
  attr(attr(dummy, "dimensions"), "raster")$affine = c(aff, 0.0)
  return(dummy)
}

s <- make_dummy_stars(6, 7, 2.5, -.5714286, -1.14) # mainly used, perspective
f <- make_dummy_stars(6, 7, 1, 1, 0) # flat
highres <- make_dummy_stars(12, 14, 1.25, -.2857143, -.57) # for resampling
lowres <- make_dummy_stars(3, 4, 5, -1, -2) # for resampling

# matrices from image ----
make_matrix <- function(image, band, n = 42, ncol = 7, t = 0) {
  # this is based on an input image with >= 4 input bands
  # n is meant to cut off NAs, ncol is y, t is random matrix for time difference
  return(matrix(image[,,,band][[1]][1:n], ncol = ncol) - t)
  # before: b3 <- matrix(warped[,,,1][[1]][1:42], ncol = 7) - t2
}

# now use function: 
b1 <- make_matrix(warped, 1)
b2 <- make_matrix(warped, 1, t = t1)
b3 <- make_matrix(warped, 1, t = t2)
g1 <- make_matrix(warped, 2)
g2 <- make_matrix(warped, 2, t = t1)
g3 <- make_matrix(warped, 2, t = t2)
r1 <- make_matrix(warped, 3)
r2 <- make_matrix(warped, 3, t = t1)
r3 <- make_matrix(warped, 3, t = t2)
n1 <- make_matrix(warped, 4)
n2 <- make_matrix(warped, 4, t = t1)
n3 <- make_matrix(warped, 4, t = t2)

# plot functions ----
plt = function(x, yoffset = 0, add, li = TRUE, pal, print_geom = TRUE, border = .75, breaks = "equal") {
  # pal is color palette
  attr(x, "dimensions")[[2]]$offset = attr(x, "dimensions")[[2]]$offset + yoffset 
  l = st_as_sf(x, as_points = FALSE)
  if (li)
    pal = lighten(pal, 0.2) # + rnorm(1, 0, 0.1))
  if (! add)
    plot(l, axes = FALSE, breaks = breaks, pal = pal, reset = FALSE, border = grey(border), key.pos = NULL, main = NULL, xlab = "time")
  else
    plot(l, axes = TRUE, breaks = breaks, pal = pal, add = TRUE, border = grey(border))
  u = st_union(l)
  # print(u)
  if(print_geom) {
    plot(st_geometry(u), add = TRUE, col = NA, border = 'black', lwd = 2.5)
  } else {
    # not print geometry
  }
}

pl_stack = function(s, x, y, add = TRUE, nrM, imgY = 7, inner = 1) {
  # nrM is the timestep {1, 2, 3}, cause this function
  # prints all 4 bands at once
  attr(s, "dimensions")[[1]]$offset = x
  attr(s, "dimensions")[[2]]$offset = y
  # m = r[[1]][y + 1:nrow,x + 1:ncol,1]
  m <- eval(parse(text=paste0("n", nrM)))
  s[[1]] = m[,c(imgY:1)] # turn around to have same orientation as flat plot
  plt(s, 0, TRUE,  pal = purples)
  m <- eval(parse(text=paste0("r", nrM)))
  s[[1]] = m[,c(imgY:1)]
  plt(s, 1*inner, TRUE,  pal = reds)
  m <- eval(parse(text=paste0("g", nrM)))
  s[[1]] = m[,c(imgY:1)]
  plt(s, 2*inner, TRUE,  pal = greens)
  m <- eval(parse(text=paste0("b", nrM)))
  s[[1]] = m[,c(imgY:1)]
  plt(s, 3*inner, TRUE, pal = blues) # li FALSE deleted
}

# flat plot function
# prints any dummy stars with any single matrix to position
pl = function(s, x, y, add = TRUE, randomize = FALSE, pal, m, print_geom = TRUE, border = .75, breaks = "equal") {
  # m is matrix to replace image with
  # m <- t(m)
  attr(s, "dimensions")[[1]]$offset = x
  attr(s, "dimensions")[[2]]$offset = y
  # print(m)
  s[[1]] = m
  plt(s, 0, add = TRUE, pal = pal, print_geom = print_geom, border = border, breaks = breaks)
  #plot(s, text_values = TRUE)
}

print_segments <- function(x, y, seg, by = 1, lwd = 4, col = "black") {
  seg = seg * by
  seg[,1] <- seg[,1] + x
  seg[,3] <- seg[,3] + x
  seg[,2] <- seg[,2] + y
  seg[,4] <- seg[,4] + y
  segments(seg[,1], seg[,2], seg[,3], seg[,4], lwd = lwd, col = col)
}

# time series ----

# from: cube1_ts_6x7_bigger.png
offset = 26
plot.new()
#par(mar = c(3, 2, 7, 2))
par(mar = c(0, 0, 0, 0))
#plot.window(xlim = c(10, 50), ylim = c(-3, 10), asp = 1)
plot.window(xlim = c(-15, 75), ylim = c(-3, 10), asp = 1)
pl_stack(s, 0, 0, nrM = 3)
pl_stack(s, offset, 0, nrM = 2)
pl_stack(s, 2 * offset, 0, nrM = 1)
# po <- matrix(c(0,-8,7,0,15,3.5,  0,1,1,5,5,14), ncol = 2)
heads <- matrix(c(3.5, 3.5 + offset, 3.5 + 2*offset, 14,14,14), ncol = 2)
points(heads, pch = 16) # 4 or 16
segments(c(-8, 7, 0, 15), c(-1,-1,3,3), 3.5, 14) # first stack pyramid
segments(c(-8, 7, 0, 15) + offset, c(-1,-1,3,3), 3.5 + offset, 14) # second stack pyramid
segments(c(-8, 7, 0, 15) + 2*offset, c(-1,-1,3,3), 3.5 + 2*offset, 14) # third stack pyramid
arrows(-13, 14, 72, 14, angle = 20, lwd = 2)  # timeline
text(7.5, 3.8, "x", col = "black")
text(-10, -2.5, "bands", srt = 90, col = "black")
text(-4.5, 1.8, "y", srt = 27.5, col = "black")
y = 15.8
text(69, y, "time", col = "black")
text(3.5, y, "2020-10-01", col = "black")
text(3.5 + offset, y, "2020-10-13", col = "black")
text(3.5 + 2*offset, y, "2020-10-25", col = "black")
# flat ----
xlabels <- seq(attr(warped, "dimensions")[[1]]$offset + attr(warped, "dimensions")[[1]]$delta / 2, length.out = attr(warped, "dimensions")[[1]]$to, by = attr(warped, "dimensions")[[1]]$delta)
ylabels <- seq(attr(warped, "dimensions")[[2]]$offset + attr(warped, "dimensions")[[2]]$delta / 2, length.out = attr(warped, "dimensions")[[2]]$to, by = attr(warped, "dimensions")[[2]]$delta)

print_labels <- function(x, y, off, lab, horizontal, cex = 1) {
  if(horizontal) { # x
    for(i in 0:(length(lab)-1)) {
      text(x + i*off, y, lab[i+1], cex = cex, srt = 90)
    }
  } else { # y
    lab <- lab[length(lab):0]
    for(i in 0:(length(lab)-1)) {
      text(x, y + i*off, lab[i+1], cex = cex)
    }
  }
}

# before: width=1000, xlim(-2, 33), date labels x=31
plot.new()
# par(mar = c(0,0,0,0))
par(mar = c(3,0,0,0))
plot.window(xlim = c(-2, 40), ylim = c(0, 25), asp = 1)
pl(f, 7, 0, pal = blues, m = b1)
pl(f, 7, 10, pal = blues, m = b2)
pl(f, 7, 20, pal = blues, m = b3)
pl(f, 14, 0, pal = greens, m = g1)
pl(f, 14, 10, pal = greens, m = g2)
pl(f, 14, 20, pal = greens, m = g3)
pl(f, 21, 0, pal = reds, m = r1)
pl(f, 21, 10, pal = reds, m = r2)
pl(f, 21, 20, pal = reds, m = r3)
pl(f, 28, 0, pal = purples, m = n1)
pl(f, 28, 10, pal = purples, m = n2)
pl(f, 28, 20, pal = purples, m = n3)
print_labels(28.5, -2, 1, xlabels, horizontal = TRUE, cex = 0.7)
print_labels(36, 0.5, 1, ylabels, horizontal = FALSE, cex = 0.7)
# arrows(6, 27, 6, 0, angle = 20, lwd = 2)
# text(5, 14, "time", srt = 90, col = "black")
text(10, 28, "blue", col = "black")
text(17, 28, "green", col = "black")
text(24, 28, "red", col = "black")
text(31, 28, "nir", col = "black")
text(3, 23.5, "2020-10-01", col = "black")
text(3, 13.5, "2020-10-13", col = "black")
text(3, 3.5, "2020-10-25", col = "black")
# filter ----
# mask <- matrix(c(rep(NA, 26), 1,NA,1,NA,1,1,1, rep(NA, 9)), ncol = 7)
mask <- matrix(c(NA,NA,NA,NA,NA,NA,
                 NA,NA,NA,NA,NA,NA,
                 NA,NA,NA, 1, 1, 1,
                 NA,NA, 1, 1, 1,NA,
                 NA,NA,NA, 1, 1,NA,
                 NA,NA,NA,NA,NA,NA,
                 NA,NA,NA,NA,NA,NA), ncol = 7)

print_grid <- function(x, y) {
  pl(f, 0+x, 0+y, pal = blues, m = b1)
  pl(f, 0+x, 10+y, pal = blues, m = b2)
  pl(f, 0+x, 20+y, pal = blues, m = b3)
  pl(f, 7+x, 0+y, pal = greens, m = g1)
  pl(f, 7+x, 10+y, pal = greens, m = g2)
  pl(f, 7+x, 20+y, pal = greens, m = g3)
  pl(f, 14+x, 0+y, pal = reds, m = r1)
  pl(f, 14+x, 10+y, pal = reds, m = r2)
  pl(f, 14+x, 20+y, pal = reds, m = r3)
  pl(f, 21+x, 0+y, pal = purples, m = n1)
  pl(f, 21+x, 10+y, pal = purples, m = n2)
  pl(f, 21+x, 20+y, pal = purples, m = n3)
}
print_alpha_grid <- function(x,y, alp = 0.2, geom = FALSE) {
  pl(f, 0+x, 0+y, pal = alpha(blues, alp), print_geom = geom,  m = b1, border = 1)
  pl(f, 0+x, 10+y, pal = alpha(blues, alp), print_geom = geom,  m = b2, border = 1)
  pl(f, 0+x, 20+y, pal = alpha(blues, alp), print_geom = geom,  m = b3, border = 1)
  pl(f, 7+x, 0+y, pal = alpha(greens, alp), print_geom = geom,  m = g1, border = 1)
  pl(f, 7+x, 10+y, pal = alpha(greens, alp), print_geom = geom,  m = g2, border = 1)
  pl(f, 7+x, 20+y, pal = alpha(greens, alp), print_geom = geom,  m = g3, border = 1)
  pl(f, 14+x, 0+y, pal = alpha(reds, alp), print_geom = geom,  m = r1, border = 1)
  pl(f, 14+x, 10+y, pal = alpha(reds, alp), print_geom = geom,  m = r2, border = 1)
  pl(f, 14+x, 20+y, pal = alpha(reds, alp), print_geom = geom,  m = r3, border = 1)
  pl(f, 21+x, 0+y, pal = alpha(purples, alp), print_geom = geom,  m = n1, border = 1)
  pl(f, 21+x, 10+y, pal = alpha(purples, alp), print_geom = geom,  m = n2, border = 1)
  invisible(pl(f, 21+x, 20+y, pal = alpha(purples, alp), print_geom = geom,  m = n3, border = 1))
}

print_grid_filter <- function(x, y) {
  pl(f, 0+x, 0+y, pal = blues, m = matrix(b1[mask == TRUE], ncol = 7))
  pl(f, 0+x, 10+y, pal = blues, m = matrix(b2[mask == TRUE], ncol = 7))
  pl(f, 0+x, 20+y, pal = blues, m = matrix(b3[mask == TRUE], ncol = 7))
  pl(f, 7+x, 0+y, pal = greens, m = matrix(g1[mask == TRUE], ncol = 7))
  pl(f, 7+x, 10+y, pal = greens, m = matrix(g2[mask == TRUE], ncol = 7))
  pl(f, 7+x, 20+y, pal = greens, m = matrix(g3[mask == TRUE], ncol = 7))
  pl(f, 14+x, 0+y, pal = reds, m = matrix(r1[mask == TRUE], ncol = 7))
  pl(f, 14+x, 10+y, pal = reds, m = matrix(r2[mask == TRUE], ncol = 7))
  pl(f, 14+x, 20+y, pal = reds, m = matrix(r3[mask == TRUE], ncol = 7))
  pl(f, 21+x, 0+y, pal = purples, m = matrix(n1[mask == TRUE], ncol = 7))
  pl(f, 21+x, 10+y, pal = purples, m = matrix(n2[mask == TRUE], ncol = 7))
  pl(f, 21+x, 20+y, pal = purples, m = matrix(n3[mask == TRUE], ncol = 7))
}

print_grid_time_filter <- function(x, y) { # 3x1, 28x7
  pl(f, 0+x, 10+y, pal = blues, m = b3)
  pl(f, 7+x, 10+y, pal = greens, m = g3)
  pl(f, 14+x, 10+y, pal = reds, m = r3)
  pl(f, 21+x, 10+y, pal = purples, m = n3)
}

print_grid_bands_filter <- function(x, y, pal = greys) { # 1x3 6x27
  pl(f, 0+x, 0+y, pal = pal, m = n1)
  pl(f, 0+x, 10+y, pal = pal, m = n2)
  pl(f, 0+x, 20+y, pal = pal, m = n3)
}

# build exactly like reduce
plot.new()
par(mar = c(3,3,3,3))
x = 120
y = 100
down = 0
plot.window(xlim = c(0, x), ylim = c(0, y), asp = 1)
print_grid(x/2-28/2,y-27)
print_alpha_grid((x/3-28)/2, 0-down) # alpha grid
print_grid_time_filter((x/3-28)/2, -10-down) # select 3rd
print_alpha_grid(x/3+((x/3-6)/2) -10.5, 0-down) # alpha grid
print_grid_bands_filter(x/3+((x/3-12.4)), 0-down, pal = purples)
print_alpha_grid(2*(x/3)+((x/3-28)/2), 0-down) # alpha grid
print_grid_filter(2*(x/3)+((x/3-28)/2), 0-down)
text(3, 13.5-down, "time", srt = 90, col = "black")
text(43, 13.5-down, "time", srt = 90, col = "black")
text(83, 13.5-down, "time", srt = 90, col = "black")
text(20, 30, "bands", col = "black")
text(60, 30, "bands", col = "black")
text(100, 30, "bands", col = "black")
arrows(x/2-28/2,y-30, x/6,32, angle = 20, lwd = 2)
arrows(x/2,y-30, x/2,32, angle = 20, lwd = 2)
arrows(x/2+28/2,y-30, 100, 32, angle = 20, lwd = 2)
# points(seq(1,120,10), seq(1,120,10))
text(28.5,49, "filter temporally", srt = 55.5, col = "black", cex = 0.8)
text(57,49, "filter bands", srt = 90, col = "black", cex = 0.8)
text(91.5,49, "filter spatially", srt = -55.5, col = "black", cex = 0.8)
print_labels(x = x/2-28/2 + 3, y = y+4, off = 7, lab = c("blue", "green", "red", "nir"),
             horizontal = TRUE, cex = 0.6)
print_labels(x = x/2-28/2 - 9, y = y-23, off = 10, lab = c("2020-10-01", "2020-10-13", "2020-10-25"),
             horizontal = FALSE, cex = 0.6)
print_labels(x = x/2-28/2 + 21.5, y = y-30, off = 1, lab = xlabels,
             horizontal = TRUE, cex = 0.3)
print_labels(x = x/2-28/2 + 30, y = y-26.5, off = 1, lab = ylabels,
             horizontal = FALSE, cex = 0.3)
# apply ----
print_text = function(s, x, y, m) {
  # m <- t(m) # transverse for correct order
  # print(m)
  r <- rep(seq(0.5, 5.5, 1), 7)
  r <- r + x # consider offsets
  u <- c(rep(0.5, 6), rep(1.5, 6), rep(2.5, 6), 
         rep(3.5, 6), rep(4.5, 6), rep(5.5, 6), rep(6.5, 6))
  u <- u + y # offset
  tab <- matrix(c(r,u), nrow = 42, byrow = FALSE) # make point table
  for (i in 1:42) {
    #text(tab[i, 1], tab[i, 2], labels = paste0("", m[i]), cex = 1.1)
    text(tab[i, 1], tab[i, 2], labels = paste0("", m[i]), cex = 0.8)
    }
}

abs_brks <- seq(-500,500, 50)
abs_pal <- sequential_hcl(n = 20, h1 = 211, c1 = 80, l1 = 30, l2 = 100, p1 = 1.2)

#png("exp_apply_unary.png", width = 2400, height = 1000, pointsize = 24)
#plot.new()
#par(mar = c(2,2,2,2))
#x = 30
#y = 7.5
#plot.window(xlim = c(0, x), ylim = c(0, y), asp = 1)
#pl(f, 3, 2, pal = abs_pal, m = b3 - 200, breaks = abs_brks)
#pl(f, 1.5, .5, pal = abs_pal, m = b2 - 200, breaks = abs_brks)
#pl(f, 0, -1, pal = abs_pal, m = b1 - 200, breaks = abs_brks)
#print_text(s, 0, -1, m = b1 - 200)
#pl(f, 23, 2, pal = abs_pal, m = abs(b3 - 200), breaks = abs_brks)
#pl(f, 21.5, 0.5, pal = abs_pal, m = abs(b2 - 200), breaks = abs_brks)
#pl(f, 20, -1, pal = abs_pal, m = abs(b1 - 200), breaks = abs_brks)
#print_text(s, 20, -1, m = abs(b1 - 200))
#arrows(11, 4, 17.5, 4, lwd = 3)
#text(14.3, 3.5, "absolute()", cex = 1.4)
#dev.off()

vNeumann_seg <- matrix(c(c(0,0,1,1,2,2,1,1,0,0,-1,-1), c(0,-1,-1,0,0,1,1,2,2,1,1,0), 
                         c(0,1,1,2,2,1,1,0,0,-1,-1,0), c(-1,-1,0,0,1,1,2,2,1,1,0,0)), ncol = 4)

apply_filter <- function(input, pad = TRUE, padValue = 1) {
  ras <- raster::focal(raster::raster(input), w = matrix(c(0,0.2,0, 0.2,0.2,0.2, 0,0.2,0), ncol = 3), pad = pad, padValue = padValue)
  ras <- raster::as.matrix(ras)
  ras[ras == "NaN"] <- -999
  return(floor(ras))
}

brks <- seq(0,1000, 50)
plot.new()
par(mar = c(0,2,0,0))
x = 30
y = 7.5
plot.window(xlim = c(0, x), ylim = c(0, y), asp = 1)
pl(f, 3, 2, pal = blues, m = b3, breaks = brks)
pl(f, 1.5, .5, pal = blues, m = b2, breaks = brks)
pl(f, 0, -1, pal = blues, m = b1, breaks = brks)
print_text(s, 0, -1, m = b1) # print text on left first stack
print_segments(2, 3, seg = vNeumann_seg, lwd = 3)
pl(f, 23, 2, pal = blues, m = apply_filter(b3), breaks = brks)
pl(f, 21.5, 0.5, pal = blues, m = apply_filter(b2), breaks = brks)
pl(f, 20, -1, pal = blues, m = apply_filter(b1), breaks = brks)
print_text(s, 20, -1, m = apply_filter(b1)) # set pad = FALSE for -99
print_segments(22, 3, seg = vNeumann_seg, lwd = 3)
arrows(11, 4, 17.5, 4, lwd = 3)
text(14.3, 3.5, "apply_kernel()", cex = 1.4)
print_segments(13.8, 1, seg = vNeumann_seg, lwd = 3)
cex = .8
text(14.3, 1.5, "0.2", cex = cex)
text(13.3, 1.5, "0.2", cex = cex)
text(15.3, 1.5, "0.2", cex = cex)
text(14.3, 2.5, "0.2", cex = cex)
text(14.3, .5, "0.2", cex = cex)

time_arrow_seg <- matrix(c(c(-1.0, 0.3, 1.5, 2.7, 3.9, 5.1), c(-1.0, 0.3, 1.5, 2.7, 3.9, 5.1),
                           c(-0.5, 0.7, 1.9, 3.1, 4.3, 5.6), c(-0.5, 0.7, 1.9, 3.1, 4.3, 5.6)), ncol = 4)
time_arrow_flag_seg <- matrix(c(c(-1.0, 1.5, 2.7, 3.9), c(-1.0, 1.5, 2.7, 3.9),
                                c(0.7, 1.9, 3.1, 5.6), c(0.7, 1.9, 3.1, 5.6)), ncol = 4)

b11 <- b2 - t1
b12 <- b1 - t2 + t1

# png("exp_apply_ts.png", width = 2400, height = 1000, pointsize = 24)
plot.new()
par(mar = c(2,2,2,2))
x = 30
y = 10 # 7.5
plot.window(xlim = c(0, x), ylim = c(0, y), asp = 1)
pl(f, 4.8, 3.8, pal = blues, m = b3, breaks = brks)
print_text(s, 4.8, 3.8, m = b3)
pl(f, 3.6, 2.6, pal = blues, m = b11, breaks = brks)
print_text(s, 3.6, 2.6, m = b11)
pl(f, 2.4, 1.4, pal = blues, m = b12, breaks = brks)
print_text(s, 2.4, 1.4, m = b12)
pl(f, 1.2, .2, pal = blues, m = b2, breaks = brks)
print_text(s, 1.2, .2, m = b2)
pl(f, 0, -1, pal = blues, m = b1, breaks = brks)
print_text(s, 0, -1, m = b1) # print text on left first stack
pl(f, 24.8, 3.8, pal = alpha(greys, 0.1), m = matrix(rep("NA", 42), ncol = 7))
pl(f, 23.6, 2.6, pal = blues, m = (b12 + b11 + b3) / 3, breaks = brks)
print_text(s, 23.6, 2.6, m = floor((b12 + b11 + b3) / 3))
pl(f, 22.4, 1.4, pal = blues, m = (b2 + b12 + b11) / 3, breaks = brks)
print_text(s, 22.4, 1.4, m = floor((b2 + b12 + b11) / 3))
pl(f, 21.2, .2, pal = blues, m = (b1 + b2 + b12) / 3, breaks = brks)
print_text(s, 21.2, .2, m = floor((b1 + b2 + b12) / 3))
pl(f, 20, -1, pal = alpha(greys, 0.1), m = matrix(rep("NA", 42), ncol = 7))
print_segments(5.7, 1.7, seg = time_arrow_seg, col = "forestgreen")
arrows(12.5, 9, 20, 9, lwd = 2)
cex = .9
text(16.3, 8.3, "apply_dimension(dimension = 't')", cex = cex)
print_segments(9.7, 1.7, time_arrow_seg, col = "forestgreen") # draw ma explanation
text(-0.5 + 10, -0.5 + 2, "496", cex = cex)
text(.7 + 10, .7 + 2, "363", cex = cex)
text(1.9 + 10, 1.9 + 2, "658", cex = cex)
text(3.1 + 10, 3.1 + 2, "230", cex = cex)
text(4.3 + 10, 4.3 + 2, "525", cex = cex)
t_formula <- expression("t"[n]*" = (t"[n-1]*" + t"[n]*" + t"[n+1]*") / 3")
# text(13.8, 3, t_formula, srt = 45, cex = 1.2)
text(14.4, 3.6, "calculate moving average", srt = 45, cex = cex)
arrows(15, 5.7, 18, 5.7, lwd = 2)
print_segments(15.4, 1.7, seg = time_arrow_seg, col = "forestgreen") # draw ma explanation
text(-0.5 + 15.7, -0.5 + 2, "NA", cex = cex)
text(.7 + 15.7, .7 + 2, "505", cex = cex)
text(1.9 + 15.7, 1.9 + 2, "417", cex = cex)
text(3.1 + 15.7, 3.1 + 2, "471", cex = cex)
text(4.3 + 15.7, 4.3 + 2, "NA", cex = cex)
print_segments(25.7, 1.7, seg = time_arrow_seg, col = "forestgreen")
# reduce ----

# calc mean over time
timeB <- (b1 + b2 + b3) / 3
timeG <- (g1 + g2 + g3) / 3
timeR <- (r1 + r2 + r3) / 3
timeN <- (n1 + n2 + n3) / 3

print_grid_time <- function(x, y) { # 3x1, 28x7
  pl(f, 0+x, 10+y, pal = blues, m = timeB)
  pl(f, 7+x, 10+y, pal = greens, m = timeG)
  pl(f, 14+x, 10+y, pal = reds, m = timeR)
  pl(f, 21+x, 10+y, pal = purples, m = timeN)
}

# calc ndvi
ndvi1 <- (n1 - r1) / (n1 + r1)
ndvi2 <- (n2 - r2) / (n2 + r2)
ndvi3 <- (n3 - r3) / (n3 + r3)

print_grid_bands <- function(x, y, pal = greys) { # 1x3 6x27
  pl(f, 0+x, 0+y, pal = pal, m = ndvi1)
  pl(f, 0+x, 10+y, pal = pal, m = ndvi2)
  pl(f, 0+x, 20+y, pal = pal, m = ndvi3)
}

plte = function(s, x, y, add = TRUE, randomize = FALSE, pal, m) {
  attr(s, "dimensions")[[1]]$offset = x
  attr(s, "dimensions")[[2]]$offset = y
  # m = r[[1]][y + 1:nrow,x + 1:ncol,1]
  # dim(m) = c(x = nrow, y = ncol) # named dim
  # s[[1]] = m
  # me <- floor(mean(s[[1]]))
  me <- floor(mean(m))
  if (me[1] > 100) { # in case non-artificial grids with very high
    me <- m / 10     # numbers are used, make them smaller
    me <- floor(mean(me))
  }
  text(x,y,me,cex = 0.8)
}

print_grid_spat <- function(x, y) {
  x = x + 3
  y = y + 3.5
  plte(s, 0+x, 0+y, pal = blues, m = b1)
  plte(s, 0+x, 10+y, pal = blues, m = b2)
  plte(s, 0+x, 20+y, pal = blues, m = b3)
  plte(s, 7+x, 0+y, pal = greens, m = g1)
  plte(s, 7+x, 10+y, pal = greens, m = g2)
  plte(s, 7+x, 20+y, pal = greens, m = g3)
  plte(s, 14+x, 0+y, pal = reds, m = r1)
  plte(s, 14+x, 10+y, pal = reds, m = r2)
  plte(s, 14+x, 20+y, pal = reds, m = r3)
  plte(s, 21+x, 0+y, pal = purples, m = n1)
  plte(s, 21+x, 10+y, pal = purples, m = n2)
  plte(s, 21+x, 20+y, pal = purples, m = n3)
}

# png("exp_reduce.png", width = 1200, height = 1000, pointsize = 32)
plot.new()
#par(mar = c(3,3,3,3))
par(mar = c(3,0,2,0))
x = 120
y = 100
plot.window(xlim = c(0, x), ylim = c(0, y), asp = 1)
print_grid(x/2-28/2,y-27)
# print_alpha_grid((x/3-28)/2, 0) # alpha grid
print_grid_time((x/3-28)/2, 0) # off = 5.5
# print_alpha_grid(x/3+((x/3-6)/2) -10.5, 0) # alpha grid
print_grid_bands(x/3+((x/3-6)/2), 0)
print_alpha_grid(2*(x/3)+((x/3-28)/2), 0, alp = 0, geom = TRUE) # alpha grid
print_grid_spat(2*(x/3)+((x/3-28)/2), 0)
text(3, 13.5, "time", srt = 90, col = "black")
#segments(3.6, 8, 3.7, 19, col = "red", lwd=3)
segments(3.4, 8, 3.4, 19, col = "red", lwd=3)
text(43, 13.5, "time", srt = 90, col = "black")
text(83, 13.5, "time", srt = 90, col = "black")
text(20, 30, "bands", col = "black")
text(60, 30, "bands", col = "black")
segments(53,29.8, 67,29.8, col = "red", lwd=3)
text(100, 30, "bands", col = "black")
text(30, 7, "x", col = "black")
text(36, 13, "y", col = "black")
text(60, -3, "x", col = "black")
text(66, 3, "y", col = "black")
text(110, -3, "x", col = "black")
text(116, 3, "y", col = "black")
segments(108,-2.4, 112,-3.2, col = "red", lwd=3)
segments(114,3.2, 118,2.4, col = "red", lwd=3)
text(60, y+4, "bands", col = "black") # dim names on main
text(43, y-14, "time", srt = 90, col = "black")
text(x/2-28/2 + 24, y-30, "x", col = "black")
text(x/2-28/2 + 30, y-24, "y", col = "black")
arrows(x/2-28/2,y-30, x/6,32, angle = 20, lwd = 2)
arrows(x/2,y-30, x/2,32, angle = 20, lwd = 2)
arrows(x/2+28/2,y-30, 100, 32, angle = 20, lwd = 2)
# points(seq(1,120,10), seq(1,120,10))
text(28.5,49, "reduce temporally", srt = 55.5, col = "black", cex = 0.8)
text(57,49, "reduce bands", srt = 90, col = "black", cex = 0.8)
text(91.5,49, "reduce spatially", srt = -55.5, col = "black", cex = 0.8)
# aggregate ----

mask_agg <- matrix(c(NA,NA,NA,NA,NA,NA,
                 NA, 1, 1,NA,NA,NA,
                  1, 1,NA,NA, 1,NA,
                  1,NA,NA, 1, 1,NA,
                  1,NA,NA, 1, 1,NA,
                  1,NA,NA,NA, 1,NA,
                 NA,NA,NA,NA,NA,NA), ncol = 7)

pl_stack_agg = function(s, x, y, add = TRUE, nrM, imgY = 7, inner = 1) {
  # pl_stack that masks the added matrices
  # nrM is the timestep {1, 2, 3}, cause this function
  # prints all 4 bands at once
  attr(s, "dimensions")[[1]]$offset = x
  attr(s, "dimensions")[[2]]$offset = y
  # m = r[[1]][y + 1:nrow,x + 1:ncol,1]
  m <- eval(parse(text=paste0("n", nrM)))
  m <- matrix(m[mask_agg == TRUE], ncol = 7)
  s[[1]] = m[,c(imgY:1)] # turn around to have same orientation as flat plot
  plt(s, 0, TRUE,  pal = purples)
  m <- eval(parse(text=paste0("r", nrM)))
  m <- matrix(m[mask_agg == TRUE], ncol = 7)
  s[[1]] = m[,c(imgY:1)]
  plt(s, 1*inner, TRUE,  pal = reds)
  m <- eval(parse(text=paste0("g", nrM)))
  m <- matrix(m[mask_agg == TRUE], ncol = 7)
  s[[1]] = m[,c(imgY:1)]
  plt(s, 2*inner, TRUE,  pal = greens)
  m <- eval(parse(text=paste0("b", nrM)))
  m <- matrix(m[mask_agg == TRUE], ncol = 7)
  s[[1]] = m[,c(imgY:1)]
  plt(s, 3*inner, TRUE, pal = blues) # li FALSE
}

polygon_1 <- matrix(c(c(0.0, 5.1, 4.9,-2.3), c(0.0, 2.4, 3.1, 1.8),
                      c(5.1, 4.9,-2.3, 0.0), c(2.4, 3.1, 1.8, 0.0)), ncol = 4)

a <- make_dummy_stars(6, 7, 5, -1.14, -2.28)

print_vector_content <- function(x, y, cex = 0.8) {
  vec <- floor(rnorm(8, 250, 100))
  text( 0 + x,12 + y, vec[1], cex = cex)
  text( 0 + x, 8 + y, vec[2], cex = cex)
  text( 0 + x, 4 + y, vec[3], cex = cex)
  text( 0 + x, 0 + y, vec[4], cex = cex)
  text( 12 + x,12 + y, vec[5], cex = cex)
  text( 12 + x, 8 + y, vec[6], cex = cex)
  text( 12 + x, 4 + y, vec[7], cex = cex)
  text( 12 + x, 0 + y, vec[8], cex = cex)
}

print_ts <- function(off2, yoff2) {
  pl_stack(s, 0 + off2, yoff2, nrM = 3) # input 2
  pl_stack(s, off + off2, yoff2, nrM = 2)
  pl_stack(s, 2 * off + off2, yoff2, nrM = 1)
  arrows(-13 + off2, 14 + yoff2, 72 + off2, 14 + yoff2, angle = 20, lwd = 2)  # timeline
  heads <- matrix(c(3.5+off2, 3.5 + off + off2, 3.5 + 2*off + off2, 14+yoff2,14+yoff2,14+yoff2), ncol = 2)
  points(heads, pch = 16) # 4 or 16
  segments(c(-8, 7, 0, 15)+off2, c(-1,-1,3,3)+yoff2, 3.5+off2, 14+yoff2) # first stack pyramid
  segments(c(-8, 7, 0, 15) + off + off2, c(-1,-1,3,3)+yoff2, 3.5 + off + off2, 14+yoff2) # second stack pyramid
  segments(c(-8, 7, 0, 15) + 2*off + off2, c(-1,-1,3,3)+yoff2, 3.5 + 2*off + off2, 14+yoff2) # third stack pyramid
  text(7.5+off2, 4.3+yoff2, "x", col = "black", cex = secText)
  text(-9.5+off2, -2.5+yoff2, "bands", srt = 90, col = "black", cex = secText)
  text(-4.5+off2, 2+yoff2, "y", srt = 27.5, col = "black", cex = secText)
  text(69+off2, 15.5+yoff2+1, "time", col = "black")
  text(3.5+off2, 15.5+yoff2, "2020-10-01", col = "black")
  text(3.5 + off + off2, 15.5+yoff2, "2020-10-13", col = "black")
  text(3.5 + 2*off + off2, 15.5+yoff2, "2020-10-25", col = "black")
}

secText = 0.8 # secondary text size (dimension naming)
off = 26 # image stacks are always 26 apart
x = 72 # png X
y = 48 # png Y
yoff = 30
plot.new()
par(mar = c(5,3,3,3))
plot.window(xlim = c(-1, x), ylim = c(0, y), asp = 1)
print_ts(5, yoff)
col = '#ff5555'
print_segments(10.57, yoff-.43, seg = polygon_1, col = col)
x = 3
segments(c(2, 0, -1.3)+x, c(-.2, 0, 1)+yoff, c(0, -1.3, 1)+x, c(0, 1, 2)+yoff, lwd = 4, col = col)
print_segments(10.57+off, yoff-.43, seg = polygon_1, col = col)
x = 3 + off
segments(c(2, 0, -1.3)+x, c(-.2, 0, 1)+yoff, c(0, -1.3, 1)+x, c(0, 1, 2)+yoff, lwd = 4, col = col)
print_segments(10.57+2*off, yoff-.43, seg = polygon_1, col = col)
x = 3 + 2 * off
segments(c(2, 0, -1.3)+x, c(-.2, 0, 1)+yoff, c(0, -1.3, 1)+x, c(0, 1, 2)+yoff, lwd = 4, col = col)
# old 75 28 poly 86.25, 27.15
pl_stack_agg(a, 5, 5, nrM = 1, inner = 3) # print masked enlargement
print_segments(16.25, 7.15, seg = polygon_1, by = 2, col = col)
x <- 1
y <- 8
segments(c(4, 0, -2.6)+x, c(-.4, 0, 2)+y, c(0, -2.6, 2)+x, c(0, 2, 4)+y, lwd = 4, col = col) # line in large
segments(-3, 25, -7, 9, lwd = 3, col = 'grey')
segments(21, 29, 27.5, 14, lwd = 3, col = 'grey')
text(10, 20, "1. Group by geometry", cex = 1.3)
vecM <- matrix(rep(1,8), ncol = 2)
text(57, 21, "2. Reduce to vector cube", cex = 1.3)
b <- make_dummy_stars(2, 4, 12, 4, 0)
pl(b, 48, -5, m = vecM, pal = alpha("white", 0.9), border = 0)
print_vector_content(54, -3)
pl(b, 46.5, -3.5, m = vecM, pal = alpha("white", 0.9), border = 0)
print_vector_content(52.5, -1.5)
pl(b, 45, -2, m = vecM, pal = alpha("white", 0.9), border = 0)
print_vector_content(51, 0)
text(51.5, 15, "Line_1", col = col)
text(63, 15, "Polygon_1", col = col)
text(57, 17.5, "Geometries", cex = 1.1)
text(42, 12, "blue")
text(42,  8, "green")
text(42,  4, "red")
text(42,  0, "nir")
text(38, 6, "Bands", srt = 90, cex = 1.1)
# arrows(13.5, -2, 13.5, -6, angle = 20, lwd = 3)
text(72, 15.5, "time", srt = 315, cex = 1.1)
arrows(69.5, 15, 72.5, 12, angle = 20, lwd = 2)
# print_segments(30, 35, seg = arrow_seg)

set.seed(1331)
library(stars)
library(colorspace)
tif = system.file("tif/L7_ETMs.tif", package = "stars")
r = read_stars(tif)

nrow = 5
ncol = 8
#m = matrix(runif(nrow * ncol), nrow = nrow, ncol = ncol)
m = r[[1]][1:nrow,1:ncol,1]
dim(m) = c(x = nrow, y = ncol) # named dim
s = st_as_stars(m)
# s
attr(s, "dimensions")[[1]]$delta = 3 
attr(s, "dimensions")[[2]]$delta = -.5
attr(attr(s, "dimensions"), "raster")$affine = c(-1.2, 0.0)

plt = function(x, yoffset = 0, add, li = TRUE) {
    attr(x, "dimensions")[[2]]$offset = attr(x, "dimensions")[[2]]$offset + yoffset 
    l = st_as_sf(x, as_points = FALSE)
    pal = sf.colors(10)
    if (li)
        pal = lighten(pal, 0.3 + rnorm(1, 0, 0.1))
    if (! add)
        plot(l, axes = FALSE, breaks = "equal", pal = pal, reset = FALSE, border = grey(.75), key.pos = NULL, main = NULL, xlab = "time")
    else
        plot(l, axes = TRUE, breaks = "equal", pal = pal, add = TRUE, border = grey(.75))
    u = st_union(l)
    plot(st_geometry(u), add = TRUE, col = NA, border = 'black', lwd = 2.5)
}

pl = function(s, x, y, add = TRUE, randomize = FALSE) {
  attr(s, "dimensions")[[1]]$offset = x
  attr(s, "dimensions")[[2]]$offset = y
  m = r[[1]][y + 1:nrow,x + 1:ncol,1]
  if (randomize)
    m = m[sample(y + 1:nrow),x + 1:ncol]
  dim(m) = c(x = nrow, y = ncol) # named dim
  s[[1]] = m
  plt(s, 0, add)
  plt(s, 1, TRUE)
  plt(s, 2, TRUE)
  plt(s, 3, TRUE)
  plt(s, 4, TRUE)
  plt(s, 5, TRUE)
  plt(s, 6, TRUE)
  plt(s, 7, TRUE)
  plt(s, 8, TRUE, FALSE)
}

# point vector data cube:
plot.new()
par(mar = c(5, 0, 5, 0))
plot.window(xlim = c(-10, 16), ylim = c(-2,12), asp = 1)
library(spacetime)
data(air)
de = st_geometry(st_normalize(st_as_sf(DE)))
# 
pl(s, 0, 0, TRUE, randomize = TRUE)
de = de * 6 + c(-7, 9)
plot(de, add = TRUE, border = grey(.5))
text(-10, 0, "time", srt = -90, col = 'black')
text(-5,  7.5, "sensor location", srt = 25, col = 'black')
text( 7,  10.5, "air quality parameter", srt = 0, col = 'black')
text( 1.5,  8.5, expression(PM[10]), col = 'black', cex = .75)
text( 4.5,  8.5, expression(NO[x]), col = 'black', cex = .75)
text( 8,  8.5, expression(SO[4]), col = 'black', cex = .75)
text( 11,  8.5, expression(O[3]), col = 'black', cex = .75)
text( 14,  8.5, expression(CO), col = 'black', cex = .75)
# location points:
p = st_coordinates(s[,1])
p[,1] = p[,1]-1.4
p[,2] = p[,2] + 8.2
points(p, col = grey(.7), pch = 16)
# centroids:
set.seed(131)
cent = st_coordinates(st_sample(de, 8))
points(cent, col = grey(.7), pch = 16)
cent = cent[rev(order(cent[,1])),]
seg = cbind(p, cent[1:8,])
segments(seg[,1], seg[,2], seg[,3], seg[,4], col = 'grey')
set.seed(1014)
options(digits = 3)

knitr::opts_chunk$set(
  comment = "#",
  collapse = knitr::is_latex_output(),
  cache = TRUE,
#  out.width = "70%",
  fig.align = 'center'
#  fig.width = 6,
#  fig.asp = 0.618,  # 1 / phi
#  fig.show = "hold"
)

options(dplyr.print_min = 6, dplyr.print_max = 6)
options(stars.crs = 17)
mapview::mapviewOptions(fgb = FALSE)

# for units: (?)
Sys.setenv(UDUNITS2_XML_PATH="")
if (knitr::is_latex_output())
    options(width = 66)
    #options(width = 72)
library(sf)
p1 = st_point(c(7.35, 52.42))
p2 = st_point(c(7.22, 52.18))
p3 = st_point(c(7.44, 52.19))
sfc = st_sfc(list(p1, p2, p3), crs = 'OGC:CRS84')
st_sf(elev = c(33.2, 52.1, 81.2), marker = c("Id01", "Id02", "Id03"),
      geom = sfc)
knitr::include_graphics("images/sf_obj.png")
library(sf)
(file = system.file("gpkg/nc.gpkg", package="sf"))
nc = st_read(file)
st_layers(file)
(file = tempfile(fileext = ".gpkg"))
st_write(nc, file, layer = "layer_nc")
nc[2:5, 3:7]
nc5 = nc[1:5, ]
nc7 = nc[1:7, ]
(i = st_intersects(nc5, nc7))
plot(st_geometry(nc7))
plot(st_geometry(nc5), add = TRUE, border = "brown")
cc = st_coordinates(st_centroid(nc7))
text(cc, labels=1:nrow(nc7), col="blue")
as.matrix(i)
lengths(i)
lengths(t(i))
methods(class = "sgbp")
library(tidyverse)
nc %>% as_tibble() %>% select(BIR74) %>% head(3)
sf_use_s2(TRUE)
orange <- nc %>% filter(NAME == "Orange")
wd = st_is_within_distance(nc, orange, units::set_units(50, km))
o50 <- nc %>% filter(lengths(wd) > 0)
nrow(o50)
og = st_geometry(orange)
plot(st_geometry(o50), lwd = 2)
plot(og, col = 'orange', add = TRUE)
plot(st_buffer(og, units::set_units(50, km)), add = TRUE, col = NA, border = 'brown')
plot(st_geometry(nc), add = TRUE, border = 'grey')
# example of largest = TRUE:
nc <- st_transform(read_sf(system.file("shape/nc.shp", package="sf")), 2264)
gr = st_sf(
         label = apply(expand.grid(1:10, LETTERS[10:1])[,2:1], 1, paste0, collapse = " "),
         geom = st_make_grid(nc))
gr$col = sf.colors(10, categorical = TRUE, alpha = .3)
# cut, to check, NA's work out:
gr = gr[-(1:30),]
suppressWarnings(nc_j <- st_join(nc, gr, largest = TRUE))
# the two datasets:
opar = par(mfrow = c(2,1), mar = rep(0,4))
plot(st_geometry(nc_j))
plot(st_geometry(gr), add = TRUE, col = gr$col)
text(st_coordinates(st_centroid(st_geometry(gr))), labels = gr$label)
# the joined dataset:
plot(st_geometry(nc_j), border = 'black', col = nc_j$col)
text(st_coordinates(st_centroid(st_geometry(nc_j))), labels = nc_j$label, cex = .8)
plot(st_geometry(gr), border = 'green', add = TRUE)
par(opar)
tif = system.file("tif/L7_ETMs.tif", package = "stars")
library(stars)
(r = read_stars(tif))
length(r)
class(r[[1]])
dim(r[[1]])
st_dimensions(r)
st_bbox(r)
tf = tempfile(fileext = ".tif")
write_stars(r, tf)
st_drivers("raster")
r[,1:100, seq(1, 250, 5), 4] %>% dim()
r[,1:100, seq(1, 250, 5), 4, drop = TRUE] %>% dim()
library(dplyr, warn.conflicts = FALSE)
filter(r, x > 289000, x < 290000)
slice(r, band, 3)
b = st_bbox(r) %>%
    st_as_sfc() %>%
    st_centroid() %>%
    st_buffer(units::set_units(500, m))
r[b]
plot(r[b][,,,1], reset = FALSE)
plot(b, border = 'brown', lwd = 2, col = NA, add = TRUE)
r[b] %>% st_normalize() %>% st_dimensions()
r[b, crop = FALSE]
st_crop(r, b)
aperm(r, c(3, 1, 2))
(rs = split(r))
merge(rs)
st_redimension(r, c(349, 352, 3, 2))
set.seed(115517)
pts = st_bbox(r) %>% st_as_sfc() %>% st_sample(20)
(e = st_extract(r, pts))
plot(r[,,,1], reset = FALSE)
col = rep("green", 20)
col[c(8, 14, 15, 18, 19)] = "red"
st_as_sf(e) %>% st_coordinates() %>% text(labels = 1:20, col = col)
rs = split(r)
trn = st_extract(rs, pts)
trn$cls = rep("land", 20)
trn$cls[c(8, 14, 15, 18, 19)] = "water"
model = MASS::lda(cls ~ ., st_set_geometry(trn, NULL))
pr = predict(rs, model)
plot(pr[1], key.pos = 4, key.width = lcm(3.5), key.length = lcm(2))
plot(r)
par(mfrow = c(1, 2))
plot(r, rgb = c(3,2,1), reset = FALSE, main = "RGB")    # rgb
plot(r, rgb = c(4,3,2), main = "False color (NIR-R-G)") # false color
log(r)
r + 2 * log(r)
r2 = r
r2[r < 50] = NA
r2
r2[is.na(r2)] = 0
r2
st_apply(r, c("x", "y"), mean)
ndvi = function(b1, b2, b3, b4, b5, b6) (b4 - b3)/(b4 + b3)
st_apply(r, c("x", "y"), ndvi)
ndvi2 = function(x) (x[4]-x[3])/(x[4]+x[3])
as.data.frame(st_apply(r, c("band"), mean))
st_apply(r, c("band"), quantile, c(.25, .5, .75))
st_apply(r, c("x", "y"), quantile, c(.25, .5, .75))
options("rgdal_show_exportToProj4_warnings"="none") # led to:
# Warning in sp::proj4string(obj): CRS object has comment, which is
# lost in output
library(spacetime)
data(air) # this loads several datasets in .GlobalEnv
dim(air)
d = st_dimensions(station = st_as_sfc(stations), time = dates)
(aq = st_as_stars(list(PM10 = air), dimensions = d))
image(aperm(log(aq), 2:1), main = "NA pattern (white) in PM10 station time series")
plot(st_as_sf(st_apply(aq, 1, mean, na.rm = TRUE)), reset = FALSE, pch = 16,
    ylim = st_bbox(DE)[c(2,4)])
plot(DE, add=TRUE)
(a = aggregate(aq, st_as_sf(DE_NUTS1), mean, na.rm = TRUE))
library(tidyverse)
a %>% filter(time >= "2008-01-01", time < "2008-01-07") %>% 
    plot(key.pos = 4)
suppressPackageStartupMessages(library(xts))
plot(as.xts(a)[,4], main = DE_NUTS1$NAME_1[4])
library(spDataLarge)
plot(st_geometry(bristol_zones), axes = TRUE, graticule = TRUE)
plot(st_geometry(bristol_zones)[33], col = 'red', add = TRUE)
head(bristol_od)
nrow(bristol_zones)^2 # all combinations
nrow(bristol_od) # non-zero combinations
# create O-D-mode array:
bristol_tidy <- bristol_od %>% 
    select(-all) %>% 
    gather("mode", "n", -o, -d)
head(bristol_tidy)
od = bristol_tidy %>% pull("o") %>% unique()
nod = length(od)
mode = bristol_tidy %>% pull("mode") %>% unique()
nmode = length(mode)
a = array(0L,  c(nod, nod, nmode), 
    dimnames = list(o = od, d = od, mode = mode))
dim(a)
a[as.matrix(bristol_tidy[c("o", "d", "mode")])] = bristol_tidy$n
order = match(od, bristol_zones$geo_code) # it happens this equals 1:102
zones = st_geometry(bristol_zones)[order]
library(stars)
(d = st_dimensions(o = zones, d = zones, mode = mode))
(odm = st_as_stars(list(N = a), dimensions = d))
plot(odm[,,33] + 1, logz = TRUE)
d = st_apply(odm, 2, sum)
which.max(d[[1]])
st_apply(odm, 1:2, sum)
st_apply(odm, c(1,3), sum)
st_apply(odm, c(2,3), sum)
o = st_apply(odm, 1, sum)
d = st_apply(odm, 2, sum)
x = (c(o, d, along = list(od = c("origin", "destination"))))
plot(x, logz = TRUE)
library(units)
a = set_units(st_area(st_as_sf(o)), km^2)
o$sum_km = o$sum / a
d$sum_km = d$sum / a
od = c(o["sum_km"], d["sum_km"], along = list(od = c("origin", "destination")))
plot(od, logz = TRUE)
sf_use_s2(TRUE)
(file = system.file("gpkg/nc.gpkg", package="sf"))
read_sf(file) %>% 
    st_geometry() %>%
    st_as_stars() %>%
    plot()
library(dplyr)
read_sf(file) %>%
    mutate(name = as.factor(NAME)) %>%
    select(SID74, SID79, name) %>%
    st_rasterize()
read_sf(file) %>%
    st_cast("MULTILINESTRING") %>%
    select(CNTY_ID) %>%
    st_rasterize() %>%
    plot()
x = st_crs("OGC:CRS84")
x$proj4string
sf_proj_search_paths()
paste0(tail(sf_proj_search_paths(), 1), .Platform$file.sep, "proj.db")
sf_extSoftVersion()["PROJ"]
sf_proj_network()
sf_proj_network(TRUE)
list.files(sf_proj_search_paths()[1], full.names = TRUE)
(p = sf_proj_pipelines("EPSG:4326", "EPSG:22525"))
sf_proj_network(FALSE)
sf_proj_pipelines("EPSG:4326", "EPSG:22525")
names(p)
p$accuracy
p[is.na(p$accuracy),]
tif = system.file("tif/L7_ETMs.tif", package = "stars")
read_stars(tif) %>%
    st_transform(4326)
read_stars(tif) %>%
    st_warp(crs = st_crs(4326)) %>%
    st_dimensions()
r = read_stars(tif)
grd = st_bbox(r) %>%
        st_as_sfc() %>%
        st_transform(4326) %>%
        st_bbox() %>%
        st_as_stars(nx = dim(r)["x"], ny = dim(r)["y"])
st_warp(r, grd)
set.seed(1014)
options(digits = 3)

knitr::opts_chunk$set(
  comment = "#",
  collapse = knitr::is_latex_output(),
  cache = TRUE,
#  out.width = "70%",
  fig.align = 'center'
#  fig.width = 6,
#  fig.asp = 0.618,  # 1 / phi
#  fig.show = "hold"
)

options(dplyr.print_min = 6, dplyr.print_max = 6)
options(stars.crs = 17)
mapview::mapviewOptions(fgb = FALSE)

# for units: (?)
Sys.setenv(UDUNITS2_XML_PATH="")
if (knitr::is_latex_output())
    options(width = 66)
    #options(width = 72)
library(sf)
(file = system.file("gpkg/nc.gpkg", package="sf"))
bb = "POLYGON ((-81.7 36.2, -80.4 36.2, -80.4 36.5, -81.7 36.5, -81.7 36.2))"
nc.1 = st_read(file, wkt_filter = bb)
q = paste("select BIR74,SID74,geom from 'nc.gpkg' where BIR74 > 1500")
nc.2 = st_read(file, query = q)
q = paste("select BIR74,SID74,geom from 'nc.gpkg' LIMIT 10 OFFSET 50")
nc.2 = st_read(file, query = q)
has_PG <- any("PostgreSQL" %in% st_drivers()$name)
nc = read_sf(file)
write_sf(nc, "PG:dbname=postgis", "nc")
pg <- DBI::dbConnect(
    RPostgres::Postgres(),
    host = "localhost",
    dbname = "postgis")
st_read(pg, query = "select BIR74,wkb_geometry from nc limit 3")
q = "SELECT BIR74,wkb_geometry FROM nc WHERE \
  ST_Intersects(wkb_geometry, 'SRID=4267;POINT (-81.49826 36.4314)');"
st_read(pg, query = q)
library(dplyr, warn.conflicts = FALSE)
nc_db = tbl(pg, "nc")
nc_db %>% 
     filter(ST_Intersects(wkb_geometry, 'SRID=4267;POINT (-81.49826 36.4314)')) %>%
     collect()
nc_db %>% filter(ST_Area(wkb_geometry) > 0.1) %>% head(3)
download.file(
  "https://openstreetmap.org/api/0.6/map?bbox=7.595,51.969,7.598,51.970"
  "data/ms.osm", method = "auto")
o = read_sf("data/ms.osm", "lines")
p = read_sf("data/ms.osm", "multipolygons")
bb = st_bbox(c(xmin=7.595, ymin = 51.969, xmax = 7.598, ymax = 51.970),
    crs = 4326)
plot(st_as_sfc(bb), axes = TRUE, lwd = 2, lty = 2, cex.axis = .5)
plot(o[,1], lwd = 2, add = TRUE)
plot(st_geometry(p), border = NA, col = '#88888888', add = TRUE)
install.packages("starsdata", repos = "http://pebesma.staff.ifgi.de", 
    type = "source")
f = "sentinel/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.zip"
granule = system.file(file = f, package = "starsdata")
file.size(granule)
base_name = strsplit(basename(granule), ".zip")[[1]]
s2 = paste0("SENTINEL2_L1C:/vsizip/", granule, "/", base_name, 
    ".SAFE/MTD_MSIL1C.xml:10m:EPSG_32632")
(p = read_stars(s2, proxy = TRUE))
object.size(p)
p2 = p * 2
plot(p)
(ds = floor(sqrt(prod(dim(p)) / prod(dev.size("px")))))
methods(class = "stars_proxy")
set.seed(1014)
options(digits = 3)

knitr::opts_chunk$set(
  comment = "#",
  collapse = knitr::is_latex_output(),
  cache = TRUE,
#  out.width = "70%",
  fig.align = 'center'
#  fig.width = 6,
#  fig.asp = 0.618,  # 1 / phi
#  fig.show = "hold"
)

options(dplyr.print_min = 6, dplyr.print_max = 6)
options(stars.crs = 17)
mapview::mapviewOptions(fgb = FALSE)

# for units: (?)
Sys.setenv(UDUNITS2_XML_PATH="")
if (knitr::is_latex_output())
    options(width = 66)
    #options(width = 72)
library(sf)
library(rnaturalearth)
#sf_use_s2(FALSE) # FIXME:
w <- ne_countries(scale = "medium", returnclass = "sf")
suppressWarnings(st_crs(w) <- st_crs(4326))
layout(matrix(1:2, 1, 2), c(2,1))
par(mar = rep(0, 4))
plot(st_geometry(w))

# sphere:
library(s2)
g = as_s2_geography(TRUE) # Earth
co = s2_data_countries()
oc = s2_difference(g, s2_union_agg(co)) # oceans
b = s2_buffer_cells(as_s2_geography("POINT(-30 -10)"), 9800000) # visible half
i = s2_intersection(b, oc) # visible ocean
co = s2_intersection(b, co)
plot(st_transform(st_as_sfc(i), "+proj=ortho +lat_0=-10 +lon_0=-30"), col = 'lightblue')
plot(st_transform(st_as_sfc(co), "+proj=ortho +lat_0=-10 +lon_0=-30"), col = NA, add = TRUE)
library(sf)
library(rnaturalearth)
w <- ne_countries(scale = "medium", returnclass = "sf")
plot(st_geometry(w))
st_is_longlat(w)
DE = st_geometry(ne_countries(country = "germany", returnclass = "sf"))
DE.eqc = st_transform(DE, "+proj=eqc +lat_ts=51.14 +lon_0=90w")
print(mean(st_bbox(DE)[c("ymin", "ymax")]), digits = 4)
par(mfrow = c(1, 2))
plot(DE, axes = TRUE)
plot(DE.eqc, axes = TRUE)
library(classInt)
# set.seed(1) needed ?
r = rnorm(100)
(cI <- classIntervals(r))
cI$brks
sf_use_s2(TRUE)
library(sf)
nc = read_sf(system.file("gpkg/nc.gpkg", package="sf"))
plot(nc["BIR74"], reset = FALSE, key.pos = 4)
plot(st_buffer(nc[1,1], units::set_units(10, km)), col = 'NA', 
     border = 'red', lwd = 2, add = TRUE)
library(stars)
r = read_stars(system.file("tif/L7_ETMs.tif", package = "stars"))
circ = st_bbox(r) %>% st_as_sfc() %>% st_sample(5) %>% st_buffer(300)
hook = function() plot(circ, col = NA, border = 'yellow', add = TRUE)
plot(r, hook = hook, key.pos = 4)
suppressPackageStartupMessages(library(tidyverse))
nc.32119 = st_transform(nc, 32119) 
year_labels = c("SID74" = "1974 - 1978", "SID79" = "1979 - 1984")
nc.32119 %>% select(SID74, SID79) %>% gather(VAR, SID, -geom) -> nc2
ggplot() + geom_sf(data = nc2, aes(fill = SID)) + 
  facet_wrap(~VAR, ncol = 1, labeller = labeller(VAR = year_labels)) +
  scale_y_continuous(breaks = 34:36) +
  scale_fill_gradientn(colors = sf.colors(20)) +
  theme(panel.grid.major = element_line(color = "white"))
library(ggplot2)
library(stars)
r = read_stars(system.file("tif/L7_ETMs.tif", package = "stars"))
ggplot() + geom_stars(data = r) +
        facet_wrap(~band) + coord_equal() +
        theme_void() +
        scale_x_discrete(expand=c(0,0)) + 
        scale_y_discrete(expand=c(0,0)) +
        scale_fill_viridis_c()
library(tmap)
tm_shape(nc.32119) + tm_polygons(c("SID74", "SID79"))
tm_shape(nc2) + tm_polygons("SID") + tm_facets(by = "VAR")
tm_shape(r) + tm_raster("L7_ETMs.tif")
tmap_mode("view")
tmap_mode("plot")
set.seed(1014)
options(digits = 3)

knitr::opts_chunk$set(
  comment = "#",
  collapse = knitr::is_latex_output(),
  cache = TRUE,
#  out.width = "70%",
  fig.align = 'center'
#  fig.width = 6,
#  fig.asp = 0.618,  # 1 / phi
#  fig.show = "hold"
)

options(dplyr.print_min = 6, dplyr.print_max = 6)
options(stars.crs = 17)
mapview::mapviewOptions(fgb = FALSE)

# for units: (?)
Sys.setenv(UDUNITS2_XML_PATH="")
if (knitr::is_latex_output())
    options(width = 66)
    #options(width = 72)
set.seed(1014)
options(digits = 3)

knitr::opts_chunk$set(
  comment = "#",
  collapse = knitr::is_latex_output(),
  cache = TRUE,
#  out.width = "70%",
  fig.align = 'center'
#  fig.width = 6,
#  fig.asp = 0.618,  # 1 / phi
#  fig.show = "hold"
)

options(dplyr.print_min = 6, dplyr.print_max = 6)
options(stars.crs = 17)
mapview::mapviewOptions(fgb = FALSE)

# for units: (?)
Sys.setenv(UDUNITS2_XML_PATH="")
if (knitr::is_latex_output())
    options(width = 66)
    #options(width = 72)
set.seed(1014)
options(digits = 3)

knitr::opts_chunk$set(
  comment = "#",
  collapse = knitr::is_latex_output(),
  cache = TRUE,
#  out.width = "70%",
  fig.align = 'center'
#  fig.width = 6,
#  fig.asp = 0.618,  # 1 / phi
#  fig.show = "hold"
)

options(dplyr.print_min = 6, dplyr.print_max = 6)
options(stars.crs = 17)
mapview::mapviewOptions(fgb = FALSE)

# for units: (?)
Sys.setenv(UDUNITS2_XML_PATH="")
if (knitr::is_latex_output())
    options(width = 66)
    #options(width = 72)
# try(load("data/ch16.rda"))
files = list.files("aq", pattern = "*.csv", full.names = TRUE)
r = lapply(files[-1], function(f) read.csv(f))
Sys.setenv(TZ = "UTC") # make sure times are not interpreted as DST
r = lapply(r, function(f) {
        f$t = as.POSIXct(f$DatetimeBegin) 
        f[order(f$t), ] 
    }
) 
r = r[sapply(r, nrow) > 1000]
names(r) =  sapply(r, function(f) unique(f$AirQualityStationEoICode))
length(r) == length(unique(names(r)))
library(xts)
r = lapply(r, function(f) xts(f$Concentration, f$t))
aq = do.call(cbind, r)
# remove stations with more than 75% missing values:
sel = apply(aq, 2, function(x) sum(is.na(x)) < 0.75 * 365 * 24)
aqsel = aq[, sel] # stations are in columns
library(tidyverse)
read.csv("aq/AirBase_v8_stations.csv", sep = "\t", stringsAsFactors = FALSE) %>% 
    as_tibble  %>% 
    filter(country_iso_code == "DE", station_type_of_area == "rural", 
                 type_of_station == "Background") -> a2
library(sf)
library(stars)
a2.sf = st_as_sf(a2, coords = c("station_longitude_deg", "station_latitude_deg"), crs = 4326)
sel =  colnames(aqsel) %in% a2$station_european_code
aqsel = aqsel[, sel]
tb = tibble(NO2 = apply(aqsel, 2, mean, na.rm = TRUE), station_european_code = colnames(aqsel))
crs = 32632
right_join(a2.sf, tb) %>% st_transform(crs) -> no2.sf 
# load German boundaries
data(air, package = "spacetime")
de <- st_transform(st_as_sf(DE_NUTS1), crs)
ggplot() + geom_sf(data = de) +  geom_sf(data = no2.sf, mapping = aes(col = NO2))
library(gstat)
v = variogram(NO2~1, no2.sf)
plot(v, plot.numbers = TRUE)
library(gstat)
v0 = variogram(NO2~1, no2.sf, cutoff = 100000, width = 10000)
plot(v0, plot.numbers = TRUE)
v.m = fit.variogram(v, vgm(1, "Exp", 50000, 1))
plot(v, v.m, plot.numbers = TRUE)
# build a grid over Germany:
st_bbox(de) %>%
  st_as_stars(dx = 10000) %>%
  st_set_crs(crs) %>%
  st_crop(de) -> grd
grd
k = krige(NO2~1, no2.sf, grd, v.m)
ggplot() + geom_stars(data = k, aes(fill = var1.pred, x = x, y = y)) + 
    geom_sf(data = st_cast(de, "MULTILINESTRING")) + 
    geom_sf(data = no2.sf)
a = aggregate(no2.sf["NO2"], by = de, FUN = mean)
b = krige(NO2~1, no2.sf, de, v.m)
b$sample = a$NO2
b$kriging = b$var1.pred
b %>% select(sample, kriging) %>% gather(var, NO2, -geometry) -> b2
ggplot() + geom_sf(data = b2, mapping = aes(fill = NO2)) + facet_wrap(~var) +
     scale_fill_gradientn(colors = sf.colors(20))
SE = function(x) sqrt(var(x)/length(x))
a = aggregate(no2.sf["NO2"], de, SE)
b$sample = a$NO2
b$kriging = sqrt(b$var1.var)
b %>% select(sample, kriging) %>% gather(var, NO2, -geometry) -> b2
ggplot() + geom_sf(data = b2, mapping = aes(fill = NO2)) + facet_wrap(~var) +
     scale_fill_gradientn(colors = sf.colors(20))
library(viridis)
s = krige(NO2~1, no2.sf, grd, v.m, nmax = 30, nsim = 10)
g = ggplot() + coord_equal() +
    scale_fill_viridis() +
    theme_void() +
    scale_x_discrete(expand=c(0,0)) +
    scale_y_discrete(expand=c(0,0))
g + geom_stars(data = s[,,,1:6]) + facet_wrap(~sample)
v = vroom::vroom("aq/pop/Zensus_Bevoelkerung_100m-Gitter.csv")
v %>% filter(Einwohner > 0) %>% 
    select(-Gitter_ID_100m) %>%
    st_as_sf(coords = c("x_mp_100m", "y_mp_100m"), crs = 3035) %>%
    st_transform(st_crs(grd)) -> b
a = aggregate(b, st_as_sf(grd, na.rm = FALSE), sum)
v = vroom::vroom("aq/pop/Zensus_Bevoelkerung_100m-Gitter.csv")
v1 <- v[v$Einwohner > 0,-1]
rm(v); gc(full=TRUE)
v2 <- st_as_sf(v1, coords = c("x_mp_100m", "y_mp_100m"), crs = 3035)
rm(v1); gc()
b <- st_transform(v2, st_crs(grd))
rm(v2); gc()
a = aggregate(b, st_as_sf(grd, na.rm = FALSE), sum)
gc()
grd$ID = 1:prod(dim(grd)) # so we can find out later which grid cell we have
ii = st_intersects(grd["ID"], st_cast(st_union(de), "MULTILINESTRING"))
grd_sf = st_as_sf(grd["ID"], na.rm = FALSE)[lengths(ii) > 0,]
iii = st_intersection(grd_sf, st_union(de))
grd$area = st_area(grd)[[1]] + units::set_units(grd$values, m^2) # NA's
grd$area[iii$ID] = st_area(iii)
grd$pop_dens = a$Einwohner / grd$area
sum(grd$pop_dens * grd$area, na.rm = TRUE) # verify
sum(b$Einwohner)
g + geom_stars(data = grd, aes(fill = sqrt(pop_dens), x = x, y = y))
(a = aggregate(grd["pop_dens"], no2.sf, mean))
no2.sf$pop_dens = st_as_sf(a)[[1]]
summary(lm(NO2~sqrt(pop_dens), no2.sf))
plot(NO2~sqrt(pop_dens), no2.sf)
abline(lm(NO2~sqrt(pop_dens), no2.sf))
no2.sf = no2.sf[!is.na(no2.sf$pop_dens),]
vr = variogram(NO2~sqrt(pop_dens), no2.sf)
vr.m = fit.variogram(vr, vgm(1, "Exp", 50000, 1))
plot(vr, vr.m, plot.numbers = TRUE)
kr = krige(NO2~sqrt(pop_dens), no2.sf, grd["pop_dens"], vr.m)
k$kr1 = k$var1.pred
k$kr2 = kr$var1.pred
st_redimension(k[c("kr1", "kr2")], 
    along = list(what = c("kriging", "residual kriging"))) %>%
    setNames("NO2") -> km
g + geom_stars(data = km, aes(fill = NO2, x = x, y = y)) + 
    geom_sf(data = st_cast(de, "MULTILINESTRING")) + 
    geom_sf(data = no2.sf) + facet_wrap(~what)
library(gstat)
demo(cokriging)
demo(cosimulation)
aqx = aq[ , colnames(aq) %in% a2$station_european_code]
sfc = st_geometry(a2.sf)[match(colnames(aqx), a2.sf$station_european_code)]
st_as_stars(NO2 = as.matrix(aqx)) %>%
    st_set_dimensions(names = c("time", "station")) %>%
    st_set_dimensions("time", index(aqx)) %>%
    st_set_dimensions("station", sfc) -> no2.st
v.st = variogramST(NO2~1, no2.st[,1:(24*31)], tlags = 0:48, 
    cores = getOption("mc.cores", 2))
v1 = plot(v.st)
v2 = plot(v.st, map = FALSE)
print(v1, split = c(1,1,2,1), more = TRUE)
print(v2, split = c(2,1,2,1), more = FALSE)
# product-sum
prodSumModel <- vgmST("productSum",
    space=vgm(150, "Exp", 200, 0),
    time= vgm(20, "Sph",   40, 0),
    k=2)
StAni = estiStAni(v.st, c(0,200000))
(fitProdSumModel <- fit.StVariogram(v.st, prodSumModel, fit.method = 7,
    stAni = StAni, method = "L-BFGS-B",
    control = list(parscale = c(1,10,1,1,0.1,1,10)),
    lower = rep(0.0001, 7)))
plot(v.st, fitProdSumModel, wireframe=FALSE, all=TRUE, scales=list(arrows=FALSE), zlim=c(0,150))
plot(v.st, model=fitProdSumModel, wireframe=TRUE, all=TRUE, scales=list(arrows=FALSE), zlim=c(0,185))
set.seed(123)
pt = st_sample(de, 2)
t = st_get_dimension_values(no2.st, 1)
st_as_stars(list(pts = matrix(1, length(t), length(pt)))) %>%
    st_set_dimensions(names = c("time", "station")) %>%
    st_set_dimensions("time", t) %>%
    st_set_dimensions("station", pt) -> new_pt
no2.st <- st_transform(no2.st, crs)
new_ts <- krigeST(NO2~1, data = no2.st["NO2"], newdata = new_pt,
         nmax = 50, stAni = StAni, modelList = fitProdSumModel,
         progress = FALSE)
plot(as.xts(new_ts[2]))
t4 = t[(1:4 - 0.5) * (3*24*30)]

d = dim(grd)
st_as_stars(pts = array(1, c(d[1], d[2], time=length(t4)))) %>%
    st_set_dimensions("time", t4) %>%
    st_set_dimensions("x", st_get_dimension_values(grd, "x")) %>%
    st_set_dimensions("y", st_get_dimension_values(grd, "y")) %>%
    st_set_crs(crs) -> grd.st
new_int <- krigeST(NO2~1, data = no2.st["NO2"], newdata = grd.st,
         nmax = 200, stAni = StAni, modelList = fitProdSumModel,
         progress = FALSE)
names(new_int)[2] = "NO2"
g + geom_stars(data = new_int, aes(fill = NO2, x = x, y = y)) + 
    facet_wrap(~time) +
    geom_sf(data = st_cast(de, "MULTILINESTRING")) + 
    geom_sf(data = no2.sf, col = 'grey', cex = .5)
set.seed(1014)
options(digits = 3)

knitr::opts_chunk$set(
  comment = "#",
  collapse = knitr::is_latex_output(),
  cache = TRUE,
#  out.width = "70%",
  fig.align = 'center'
#  fig.width = 6,
#  fig.asp = 0.618,  # 1 / phi
#  fig.show = "hold"
)

options(dplyr.print_min = 6, dplyr.print_max = 6)
options(stars.crs = 17)
mapview::mapviewOptions(fgb = FALSE)

# for units: (?)
Sys.setenv(UDUNITS2_XML_PATH="")
if (knitr::is_latex_output())
    options(width = 66)
    #options(width = 72)
library(sf)
data(pol_pres15, package="spDataLarge")
head(pol_pres15[, c(1, 4, 6)])
library(tmap)
tm_shape(pol_pres15) + tm_fill("types")
library(spdep)
args(poly2nb)
system.time(nb_q <- poly2nb(pol_pres15, queen=TRUE))
nb_q
system.time(nb_q_legacy <- poly2nb(pol_pres15, queen=TRUE, small_n=2500))
all.equal(nb_q, nb_q_legacy, check.attributes=FALSE)
n.comp.nb(nb_q)$nc
tf <- tempfile(fileext=".gal")
write.nb.gal(nb_q, tf)
coords <- st_centroid(st_geometry(pol_pres15), of_largest_polygon=TRUE)
suppressMessages(nb_tri <- tri2nb(coords))
nb_tri
summary(unlist(nbdists(nb_tri, coords)))
n.comp.nb(nb_tri)$nc
nb_soi <- graph2nb(soi.graph(nb_tri, coords))
nb_soi
n_comp <- n.comp.nb(nb_soi)
n_comp$nc
table(n_comp$comp.id)
opar <- par(mar=c(0,0,0,0)+0.5)
plot(st_geometry(pol_pres15), border="grey", lwd=0.5)
plot(nb_soi, coords=st_coordinates(coords), add=TRUE, points=FALSE, lwd=0.5)
plot(diffnb(nb_tri, nb_soi), coords=st_coordinates(coords), col="orange",
     add=TRUE, points=FALSE, lwd=0.5)
par(opar)
k1 <- knn2nb(knearneigh(coords))
k1dists <- unlist(nbdists(k1, coords))
summary(k1dists)
system.time(nb_d18 <- dnearneigh(coords, 0, 18000))
system.time(nb_d18a <- dnearneigh(coords, 0, 18000, use_kd_tree=FALSE))
all.equal(nb_d18, nb_d18a, check.attributes=FALSE)
nb_d18
n_comp <- n.comp.nb(nb_d18)
n_comp$nc
table(n_comp$comp.id)
nb_d183 <- dnearneigh(coords, 0, 18300)
nb_d183
n_comp <- n.comp.nb(nb_d183)
n_comp$nc
nb_d16 <- dnearneigh(coords, 0, 16000)
nb_d16
knn_k6 <- knearneigh(coords, k=6)
knn2nb(knn_k6)
nb_k6s <- knn2nb(knn_k6, sym=TRUE)
nb_k6s
n_comp <- n.comp.nb(nb_k6s)
n_comp$nc
args(nb2listw)
args(spweights.constants)
lw_q_B <- nb2listw(nb_q, style="B")
unlist(spweights.constants(lw_q_B))
lw_q_W <- nb2listw(nb_q, style="W")
unlist(spweights.constants(lw_q_W))[c(1,6:8)]
gwts <- lapply(nbdists(nb_d183, coords), function(x) 1/(x/1000))
lw_d183_idw_B <- nb2listw(nb_d183, glist=gwts, style="B")
unlist(spweights.constants(lw_d183_idw_B))[c(1,6:8)]
try(lw_d16_B <- nb2listw(nb_d16, style="B"))
lw_d16_B <- nb2listw(nb_d16, style="B", zero.policy=TRUE)
unlist(spweights.constants(lw_d16_B,  zero.policy=TRUE))[c(1,6:8)]
set.seed(1)
x <- rnorm(nrow(pol_pres15))
mt <- moran.test(x, lw_q_B, randomisation=FALSE, alternative="two.sided")
glance_htest <- function(ht) c(ht$estimate, "Std deviate"=unname(ht$statistic), "p.value"=unname(ht$p.value))
glance_htest(mt)
beta <- 0.15e-02
t <- st_coordinates(coords)[,1]/1000
x_t <- x + beta*t
mt <- moran.test(x_t, lw_q_B, randomisation=FALSE, alternative="two.sided")
glance_htest(mt)
lmt <- lm.morantest(lm(x_t ~ t), lw_q_B, alternative="two.sided")
glance_htest(lmt)
args(joincount.test)
table(pol_pres15$types)
joincount.multi(pol_pres15$types, listw=lw_q_B)
joincount.multi(pol_pres15$types, listw=lw_d183_idw_B)
args(moran.test)
mt <- moran.test(pol_pres15$I_turnout, listw=lw_q_B, randomisation=FALSE)
glance_htest(mt)
ols <- lm(I_turnout ~ 1, pol_pres15)
lmt <- lm.morantest(ols, listw=lw_q_B)
glance_htest(lmt)
mtr <- moran.test(pol_pres15$I_turnout, listw=lw_q_B)
glance_htest(mtr)
set.seed(1)
mmc <- moran.mc(pol_pres15$I_turnout, listw=lw_q_B, nsim=999, return_boot = TRUE)
c("Permutation bootstrap"=var(mmc$t), 
  "Analytical randomisation"=unname(mtr$estimate[3]))
infl_W <- moran.plot(pol_pres15$I_turnout, listw=lw_q_W, labels=pol_pres15$TERYT, cex=1, pch=".", xlab="I round turnout", ylab="lagged turnout")
pol_pres15$hat_value <- infl_W$hat
tm_shape(pol_pres15) + tm_fill("hat_value")
locm <- localmoran(pol_pres15$I_turnout, listw=lw_q_W, alternative="two.sided")
all.equal(sum(locm[,1])/Szero(lw_q_W), 
          unname(moran.test(pol_pres15$I_turnout, lw_q_W)$estimate[1]))
pvs <- cbind("none"=locm[,5], "bonferroni"=p.adjust(locm[,5], "bonferroni"), 
             "fdr"=p.adjust(locm[,5],"fdr"), "BY"=p.adjust(locm[,5], "BY"))
apply(pvs, 2, function(x) sum(x < 0.05))
library(parallel)
set.coresOption(ifelse(detectCores() == 1, 1, detectCores()-1L))
system.time(locm_p <- localmoran_perm(pol_pres15$I_turnout, listw=lw_q_W, nsim=499, alternative="two.sided", iseed=1))
pv <- locm_p[,5]
pvsp <- cbind("none"=pv, "bonferroni"=p.adjust(pv, "bonferroni"), "fdr"=p.adjust(pv, "fdr"), "BY"=p.adjust(pv, "BY"))
apply(pvsp, 2, function(x) sum(x <= 0.05))
brks <- qnorm(c(0, 0.00001, 0.0001, 0.001, 0.01, 0.025, 0.5, 0.975, 0.99, 0.999, 0.9999, 0.99999, 1))
(tab <- table(cut(locm_p[,4], brks)))
sum(tab[c(1:5, 8:12)])
sum(tab[c(1, 12)])
pol_pres15$locm_Z <- locm[,4]
pol_pres15$locm_p_Z <- locm_p[,4]
tm_shape(pol_pres15) + tm_fill(c("locm_Z", "locm_p_Z"), breaks=brks, midpoint=0, title="Standard deviates of\nLocal Moran's I") + tm_facets(free.scales=FALSE) + tm_layout(panel.labels=c("Analytical", "Conditional permutation"))
quadr <- interaction(cut(infl_W$x, c(-Inf, mean(infl_W$x), Inf), labels=c("Low X", "High X")), cut(infl_W$wx, c(-Inf, mean(infl_W$wx), Inf), labels=c("Low WX", "High WX")), sep=" : ")
a <- table(quadr)
pol_pres15$hs_an_q <- quadr
is.na(pol_pres15$hs_an_q) <- !(pol_pres15$locm_Z < brks[6] | pol_pres15$locm_Z > brks[8])
b <- table(pol_pres15$hs_an_q)
pol_pres15$hs_cp_q <- quadr
is.na(pol_pres15$hs_cp_q) <- !(pol_pres15$locm_p_Z < brks[2] | pol_pres15$locm_p_Z > brks[12])
c <- table(pol_pres15$hs_cp_q)
t(rbind("Moran plot quadrants"=a, "Unadjusted analytical"=b, "Bonferroni cond. perm."=c))
tm_shape(pol_pres15) + tm_fill(c("hs_an_q", "hs_cp_q"), colorNA="grey95", textNA="Not significant", title="Turnout hotspot status\nLocal Moran's I") + tm_facets(free.scales=FALSE) + tm_layout(panel.labels=c("Unadjusted analytical", "Cond. perm. with Bonferroni"))
system.time(locG <- localG(pol_pres15$I_turnout, lw_q_W))
system.time(locG_p <- localG_perm(pol_pres15$I_turnout, lw_q_W, nsim=499, iseed=1))
pv <- 2 * pnorm(abs(c(locG)), lower.tail = FALSE)
pvsp <- cbind(pv, p.adjust(pv, "bonferroni"), p.adjust(pv, "fdr"), 
              p.adjust(pv, "BY"))
colnames(pvsp) <- c("none", "bonferroni", "fdr", "BY")
apply(pvsp, 2, function(x) sum(x < 0.05))
library(ggplot2)
p1 <- ggplot(data.frame(Zi=locm[,4], Zi_perm=locm_p[,4])) + geom_point(aes(x=Zi, y=Zi_perm), alpha=0.2) + xlab("Analytical") + ylab("Conditional permutation") + coord_fixed() + ggtitle("Local Moran's I")
p2 <- ggplot(data.frame(Zi=c(locG), Zi_perm=c(locG_p))) + geom_point(aes(x=Zi, y=Zi_perm), alpha=0.2) + xlab("Analytical") + ylab("Conditional permutation") + coord_fixed() + ggtitle("Local G")
gridExtra::grid.arrange(p1, p2, nrow=1)
pol_pres15$locG_Z <- c(locG)
pol_pres15$hs_G <- cut(c(locG), c(-Inf, brks[2], brks[12], Inf), labels=c("Low", "Not significant", "High"))
table(pol_pres15$hs_G)
m1 <- tm_shape(pol_pres15) + tm_fill(c("locG_Z"), midpoint=0, title="Standard\ndeviate")
m2 <- tm_shape(pol_pres15) + tm_fill(c("hs_G"), title="Bonferroni\nhotspot status")
tmap_arrange(m1, m2, nrow=1)
library(rgeoda)
system.time(Geoda_w <- queen_weights(pol_pres15))
summary(Geoda_w)
system.time(lisa <- local_moran(Geoda_w, pol_pres15["I_turnout"], 
    cpu_threads=ifelse(parallel::detectCores() == 1, 1, parallel::detectCores()-1L), permutations=499, seed=1))
all.equal(card(nb_q), lisa_num_nbrs(lisa), check.attributes=FALSE)
all.equal(lisa_values(lisa), localmoran(pol_pres15$I_turnout, listw=lw_q_W, mlvar=FALSE)[,1], check.attributes=FALSE)
set.seed(1014)
options(digits = 3)

knitr::opts_chunk$set(
  comment = "#",
  collapse = knitr::is_latex_output(),
  cache = TRUE,
#  out.width = "70%",
  fig.align = 'center'
#  fig.width = 6,
#  fig.asp = 0.618,  # 1 / phi
#  fig.show = "hold"
)

options(dplyr.print_min = 6, dplyr.print_max = 6)
options(stars.crs = 17)
mapview::mapviewOptions(fgb = FALSE)

# for units: (?)
Sys.setenv(UDUNITS2_XML_PATH="")
if (knitr::is_latex_output())
    options(width = 66)
    #options(width = 72)
library(sf)
library(spatialreg)
boston_506 <- st_read(system.file("shapes/boston_tracts.shp", package="spData")[1])
nb_q <- spdep::poly2nb(boston_506)
lw_q <- spdep::nb2listw(nb_q, style="W")
table(boston_506$censored)
summary(boston_506$median)
boston_489 <- boston_506[!is.na(boston_506$median),]
nb_q_489 <- spdep::poly2nb(boston_489)
lw_q_489 <- spdep::nb2listw(nb_q_489, style="W", zero.policy=TRUE)
agg_96 <- list(as.character(boston_506$NOX_ID))
boston_96 <- aggregate(boston_506[, "NOX_ID"], by=agg_96, unique)
nb_q_96 <- spdep::poly2nb(boston_96)
lw_q_96 <- spdep::nb2listw(nb_q_96)
boston_96$NOX <- aggregate(boston_506$NOX, agg_96, mean)$x
boston_96$CHAS <- aggregate(as.integer(boston_506$CHAS)-1, agg_96, max)$x
nms <- names(boston_506)
ccounts <- 23:31
for (nm in nms[c(22, ccounts, 36)]) {
  boston_96[[nm]] <- aggregate(boston_506[[nm]], agg_96, sum)$x
}
br2 <- c(3.50,  6.25,  8.75, 12.50, 17.50, 22.50, 30.00, 42.50, 60.00)*1000
counts <- as.data.frame(boston_96)[, nms[ccounts]]
f <- function(x) matrixStats::weightedMedian(x=br2, w=x, interpolate=TRUE)
boston_96$median <- apply(counts, 1, f)
is.na(boston_96$median) <- boston_96$median > 50000
summary(boston_96$median)
POP <- boston_506$POP
f <- function(x) matrixStats::weightedMean(x[,1], x[,2])
for (nm in nms[c(9:11, 14:19, 21, 33)]) {
  s0 <- split(data.frame(boston_506[[nm]], POP), agg_96)
  boston_96[[nm]] <- sapply(s0, f)
}
boston_94 <- boston_96[!is.na(boston_96$median),]
nb_q_94 <- spdep::subset.nb(nb_q_96, !is.na(boston_96$median))
lw_q_94 <- spdep::nb2listw(nb_q_94, style="W")
form <- formula(log(median) ~ CRIM + ZN + INDUS + CHAS + I((NOX*10)^2) + I(RM^2) + 
                  AGE + log(DIS) + log(RAD) + TAX + PTRATIO + I(BB/100) + 
                  log(I(LSTAT/100)))
args(errorsarlm)
args(lagsarlm)
args(spautolm)
args(sacsarlm)
args(getS3method("summary", "Sarlm"))
eigs_489 <- eigenw(lw_q_489)
SDEM_489 <- errorsarlm(form, data=boston_489, listw=lw_q_489, Durbin=TRUE, 
                       zero.policy=TRUE, control=list(pre_eig=eigs_489))
SEM_489 <- errorsarlm(form, data=boston_489, listw=lw_q_489, 
                      zero.policy=TRUE, control=list(pre_eig=eigs_489))
cbind(data.frame(model=c("SEM", "SDEM")), 
      rbind(broom::tidy(Hausman.test(SEM_489)), 
            broom::tidy(Hausman.test(SDEM_489))))[,1:4]
cbind(data.frame(model=c("SEM", "SDEM")), 
      rbind(broom::tidy(LR1.Sarlm(SEM_489)), 
            broom::tidy(LR1.Sarlm(SDEM_489))))[,c(1, 4:6)]
eigs_94 <- eigenw(lw_q_94)
SDEM_94 <- errorsarlm(form, data=boston_94, listw=lw_q_94, Durbin=TRUE,
                      control=list(pre_eig=eigs_94))
SEM_94 <- errorsarlm(form, data=boston_94, listw=lw_q_94, control=list(pre_eig=eigs_94))
cbind(data.frame(model=c("SEM", "SDEM")), 
      rbind(broom::tidy(Hausman.test(SEM_94)), 
            broom::tidy(Hausman.test(SDEM_94))))[, 1:4]
cbind(data.frame(model=c("SEM", "SDEM")), rbind(broom::tidy(LR1.Sarlm(SEM_94)), broom::tidy(LR1.Sarlm(SDEM_94))))[,c(1, 4:6)]
broom::tidy(lmtest::lrtest(SEM_489, SDEM_489))
broom::tidy(lmtest::lrtest(SEM_94, SDEM_94))
args(lmSLX)
SLX_489 <- lmSLX(form, data=boston_489, listw=lw_q_489, zero.policy=TRUE)
broom::tidy(lmtest::lrtest(SLX_489, SDEM_489))
SLX_94 <- lmSLX(form, data=boston_94, listw=lw_q_94)
broom::tidy(lmtest::lrtest(SLX_94, SDEM_94))
SLX_94w <- lmSLX(form, data=boston_94, listw=lw_q_94, weights=units)
SDEM_94w <- errorsarlm(form, data=boston_94, listw=lw_q_94, Durbin=TRUE, weights=units,
                       control=list(pre_eig=eigs_94))
broom::tidy(lmtest::lrtest(SLX_94w, SDEM_94w))
SEM1_94 <- sphet::spreg(form, data=boston_94, listw=lw_q_94, model="error")
res <- rbind(summary(SEM_94)$Coef["I((NOX * 10)^2)",], 
             summary(SEM1_94)$CoefTable["I((NOX * 10)^2)",])
rownames(res) <- c("ML", "GMM")
res
formiv <- update(form, . ~ . - I((NOX*10)^2))
boston_94$NOX2 <- (boston_94$NOX*10)^2
suppressWarnings(ccoords <- st_coordinates(st_centroid(st_geometry(boston_94))))
iform <- formula(~poly(ccoords, degree=2) + DIS + RAD)
SEM1_94iv <- sphet::spreg(formiv, data=boston_94, listw=lw_q_94, endog = ~NOX2, 
                   instruments=iform, model="error")
summary(SEM1_94iv)$CoefTable["NOX2",]
system.time(SDEM_94B <- spBreg_err(form, data=boston_94, listw=lw_q_94, Durbin=TRUE))
system.time(SDEM_489B <- spBreg_err(form, data=boston_489, listw=lw_q_489, Durbin=TRUE, zero.policy=TRUE))
t(errorsarlm(form, data=boston_94, listw=lw_q_94, Durbin=TRUE)$timings[,2])
t(errorsarlm(form, data=boston_489, listw=lw_q_489, Durbin=TRUE, zero.policy=TRUE)$timings[,2])
t(attr(SDEM_94B, "timings")[ , 3])
t(attr(SDEM_489B, "timings")[ , 3])
SEM_94$interval
1/range(eigs_94)
1/c(lextrW(lw_q_94))
W <- as(lw_q_94, "CsparseMatrix")
1/Re(c(RSpectra::eigs(W, k=1, which="SR")$values, 
       RSpectra::eigs(W, k=1, which="LR")$values))
coef <- 0.5
sum(log(1 - coef * eigs_94))
I <- Diagonal(nrow(boston_94))
LU <- lu(I - coef * W)
dU <- abs(diag(slot(LU, "U")))
sum(log(dU))
W <- as(similar.listw(lw_q_94), "CsparseMatrix")
super <- as.logical(NA)
cch <- Cholesky((I - coef * W), super=super)
c(2 * determinant(cch, logarithm = TRUE)$modulus)
args(jacobianSetup)
args(do_ldet)
args(getS3method("impacts", "Sarlm"))
args(getS3method("summary", "LagImpact"))
sum_imp_94_SDEM <- summary(impacts(SDEM_94))
rbind(Impacts=sum_imp_94_SDEM$mat[5,], SE=sum_imp_94_SDEM$semat[5,])
sum_imp_94_SDEM_B <- summary(impacts(SDEM_94B))
rbind(Impacts=sum_imp_94_SDEM_B$mat[5,], SE=sum_imp_94_SDEM_B$semat[5,])
sum_imp_94_SLX <- summary(impacts(SLX_94))
rbind(Impacts=sum_imp_94_SLX$mat[5,], SE=sum_imp_94_SLX$semat[5,])
SLM_489 <- lagsarlm(form, data=boston_489, listw=lw_q_489, zero.policy=TRUE)
args(trW)
W <- as(lw_q_489, "CsparseMatrix")
tr_489 <- trW(W)
str(tr_489)
SLM_489_imp <- impacts(SLM_489, tr=tr_489, R=2000)
SLM_489_imp_sum <- summary(SLM_489_imp, short=TRUE, zstats=TRUE)
res <- rbind(Impacts=sapply(SLM_489_imp$res, "[", 5), SE=SLM_489_imp_sum$semat[5,])
colnames(res) <- c("Direct", "Indirect", "Total")
res
coef_SLM_489 <- coef(SLM_489)
IrW <- Diagonal(489) - coef_SLM_489[1] * W
S_W <- solve(IrW)
S_NOX_W <- S_W %*% (diag(489) * coef_SLM_489[7])
c(Direct=mean(diag(S_NOX_W)), Total=sum(S_NOX_W)/489)
sapply(impacts(SLM_489, listw=lw_q_489), "[", 5)
sapply(impacts(SLM_489, evalues=eigs_489), "[", 5)
args(getS3method("predict", "sarlm"))
nd_489 <- boston_489
nd_489$PTRATIO <- nd_489$PTRATIO + 1
OLS_489 <- lm(form, data=boston_489)
fitted <- predict(OLS_489)
nd_fitted <- predict(OLS_489, newdata=nd_489)
all.equal(unname(coef(OLS_489)[12]), mean(nd_fitted - fitted))
fitted <- predict(SLM_489)
nd_fitted <- predict(SLM_489, newdata=nd_489, listw=lw_q_489, pred.type="TS",
                     zero.policy=TRUE)
all.equal(unname(coef_SLM_489[13]), mean(nd_fitted - fitted))
nd <- boston_506[is.na(boston_506$median),]
t0 <- exp(predict(SDEM_489, newdata=nd, listw=lw_q, pred.type="TS", zero.policy=TRUE))
suppressWarnings(t1  <- exp(predict(SDEM_489, newdata=nd, listw=lw_q, pred.type="KP2",
                                    zero.policy=TRUE)))
suppressWarnings(t2  <- exp(predict(SDEM_489, newdata=nd, listw=lw_q, pred.type="KP5",
                                    zero.policy=TRUE)))
data.frame(fit_TS=t0[,1], fit_KP2=c(t1), fit_KP5=c(t2),
           censored=boston_506$censored[as.integer(attr(t0, "region.id"))])
boston_94a <- aggregate(boston_489[,"NOX_ID"], list(boston_489$NOX_ID), unique)
nb_q_94a <- spdep::poly2nb(boston_94a)
NOX_ID_no_neighs <- boston_94a$NOX_ID[which(spdep::card(nb_q_94a) == 0)]
boston_487 <- boston_489[is.na(match(boston_489$NOX_ID, NOX_ID_no_neighs)),]
boston_93 <- aggregate(boston_487[, "NOX_ID"], list(ids = boston_487$NOX_ID), unique)
row.names(boston_93) <- as.character(boston_93$NOX_ID)
nb_q_93 <- spdep::poly2nb(boston_93, row.names=unique(as.character(boston_93$NOX_ID)))
library(lme4)
MLM <- lmer(update(form, . ~ . + (1 | NOX_ID)), data=boston_487, REML=FALSE)
boston_93$MLM_re <- ranef(MLM)[[1]][,1]
library(Matrix)
suppressMessages(library(MatrixModels))
Delta <- as(model.Matrix(~ -1 + as.factor(NOX_ID), data=boston_487, sparse=TRUE),
            "dgCMatrix")
M <- as(spdep::nb2listw(nb_q_93, style="B"), "CsparseMatrix")
suppressPackageStartupMessages(library(hglm))
y_hglm <- log(boston_487$median)
X_hglm <- model.matrix(lm(form, data=boston_487))
suppressWarnings(HGLM_iid <- hglm(y=y_hglm, X=X_hglm, Z=Delta))
suppressWarnings(HGLM_sar <- hglm(y=y_hglm, X=X_hglm, Z=Delta, rand.family=SAR(D=M)))
boston_93$HGLM_re <- unname(HGLM_iid$ranef)
boston_93$HGLM_ss <- HGLM_sar$ranef[,1]
library(HSAR)
suppressWarnings(HSAR <- hsar(form, data=boston_487, W=NULL, M=M, Delta=Delta, 
                              burnin=500, Nsim=2500, thinning=1))
boston_93$HSAR_ss <- HSAR$Mus[1,]
suppressPackageStartupMessages(library(R2BayesX))
BX_iid <- bayesx(update(form, . ~ . + sx(NOX_ID, bs="re")), family="gaussian",
data=boston_487, method="MCMC", iterations=12000, burnin=2000, step=2, seed=123)
boston_93$BX_re <- BX_iid$effects["sx(NOX_ID):re"][[1]]$Mean
RBX_gra <- nb2gra(nb_q_93)
BX_mrf <- bayesx(update(form, . ~ . + sx(NOX_ID, bs="mrf", map=RBX_gra)), 
                 family="gaussian", data=boston_487, method="MCMC", 
                 iterations=12000, burnin=2000,step=2, seed=123)
boston_93$BX_ss <- BX_mrf$effects["sx(NOX_ID):mrf"][[1]]$Mean
library(mgcv)
names(nb_q_93) <- attr(nb_q_93, "region.id")
boston_487$NOX_ID <- as.factor(boston_487$NOX_ID)
GAM_MRF <- gam(update(form, . ~ . + s(NOX_ID, bs="mrf", xt=list(nb=nb_q_93))),
               data=boston_487, method="REML")
boston_93$GAM_ss <- aggregate(predict(GAM_MRF, type="terms", se=FALSE)[,14],
                              list(boston_487$NOX_ID), mean)$x
res <- rbind(iid_lmer=summary(MLM)$coefficients[6, 1:2],
             iid_hglm=summary(HGLM_iid)$FixCoefMat[6, 1:2], 
             iid_BX=BX_iid$fixed.effects[6, 1:2], 
             sar_hsar=c(HSAR$Mbetas[1, 6], HSAR$SDbetas[1, 6]),
             mrf_BX=BX_mrf$fixed.effects[6, 1:2], 
             mrf_GAM=c(summary(GAM_MRF)$p.coeff[6], summary(GAM_MRF)$se[6]))
suppressPackageStartupMessages(library(ggplot2))
df_res <- as.data.frame(res)
names(df_res) <- c("mean", "sd")
limits <- aes(ymax = mean + qnorm(0.975)*sd, ymin=mean + qnorm(0.025)*sd)
df_res$model <- row.names(df_res)
p <- ggplot(df_res, aes(y=mean, x=model)) + geom_point() + geom_errorbar(limits) + geom_hline(yintercept = 0, col="#EB811B") + coord_flip()
p + ggtitle("NOX coefficients and error bars") + theme(plot.background = element_rect(fill = "transparent",colour = NA), legend.background = element_rect(colour = NA, fill = "transparent"))
library(tmap)
# FIXME: uncommented BX_re
tm_shape(boston_93) + tm_fill(c("MLM_re", "HGLM_re" #, "BX_re"
                                ), midpoint=0, title="IID")  + tm_facets(free.scales=FALSE) + tm_borders(lwd=0.3, alpha=0.4) + tm_layout(panel.labels=c("MLM", "HGLM", "BX"))
# FIXME: uncommented BX_ss
tm_shape(boston_93) + tm_fill(c("HGLM_ss", "HSAR_ss", # "BX_ss", 
  "GAM_ss"), midpoint=0, title="SS")  + tm_facets(free.scales=FALSE) + tm_borders(lwd=0.3, alpha=0.4) + tm_layout(panel.labels=c("HGLM SAR", "HSAR SAR", "BX MRF", "GAM MRF"))
set.seed(1014)
options(digits = 3)

knitr::opts_chunk$set(
  comment = "#",
  collapse = knitr::is_latex_output(),
  cache = TRUE,
#  out.width = "70%",
  fig.align = 'center'
#  fig.width = 6,
#  fig.asp = 0.618,  # 1 / phi
#  fig.show = "hold"
)

options(dplyr.print_min = 6, dplyr.print_max = 6)
options(stars.crs = 17)
mapview::mapviewOptions(fgb = FALSE)

# for units: (?)
Sys.setenv(UDUNITS2_XML_PATH="")
if (knitr::is_latex_output())
    options(width = 66)
    #options(width = 72)
set.seed(1014)
options(digits = 3)

knitr::opts_chunk$set(
  comment = "#",
  collapse = knitr::is_latex_output(),
  cache = TRUE,
#  out.width = "70%",
  fig.align = 'center'
#  fig.width = 6,
#  fig.asp = 0.618,  # 1 / phi
#  fig.show = "hold"
)

options(dplyr.print_min = 6, dplyr.print_max = 6)
options(stars.crs = 17)
mapview::mapviewOptions(fgb = FALSE)

# for units: (?)
Sys.setenv(UDUNITS2_XML_PATH="")
if (knitr::is_latex_output())
    options(width = 66)
    #options(width = 72)
# All R code in this book {-}
set.seed(1014)
options(digits = 3)

knitr::opts_chunk$set(
  comment = "#",
  collapse = knitr::is_latex_output(),
  cache = TRUE,
#  out.width = "70%",
  fig.align = 'center'
#  fig.width = 6,
#  fig.asp = 0.618,  # 1 / phi
#  fig.show = "hold"
)

options(dplyr.print_min = 6, dplyr.print_max = 6)
options(stars.crs = 17)
mapview::mapviewOptions(fgb = FALSE)

# for units: (?)
Sys.setenv(UDUNITS2_XML_PATH="")
if (knitr::is_latex_output())
    options(width = 66)
    #options(width = 72)
a %>% b() %>% c() %>% d(n = 10)
d(c(b(a)), n = 10)
tmp1 <- b(a)
tmp2 <- c(tmp1)
tmp3 <- d(tmp2, n = 10)
typeof(1:10)
length(1:10)
typeof(1.0)
length(1.0)
typeof(c("foo", "bar"))
length(c("foo", "bar"))
typeof(c(TRUE, FALSE))
i = integer(0)
typeof(i)
i
length(i)
a = c(1,2,3)
a[2]
a[[2]]
a[2:3]
a[2:3] = c(5,6)
a
a[[3]] = 10
a
l <- list(3, TRUE, "foo")
typeof(l)
length(l)
l[1]
l[[1]]
l[1:2] = list(4, FALSE)
l
l[[3]] = "bar"
l
l = list(first = 3, second = TRUE, third = "foo")
l
l$second
l$second = FALSE
l
3 == NULL # not FALSE!
NULL == NULL # not even TRUE!
is.null(NULL)
l = l[c(1,3)] # remove second, implicitly
l
l$second = NULL
l
a = 1:3
attr(a, "some_meta_data") = "foo"
a
attr(a, "some_meta_data")
attr(a, "some_meta_data") = "bar"
attr(a, "some_meta_data")
attributes(a)
attributes(a) = list(some_meta_data = "foo")
attributes(a)
class(1:3)
class(c(TRUE, FALSE))
class(c("TRUE", "FALSE"))
a = 1:3
class(a) = "foo"
a
class(a)
attributes(a)
print.foo = function(x, ...) print(paste("an object of class foo with length", length(x)))
print(a)
unclass(a)
library(sf)
p = st_polygon(list(rbind(c(0,0), c(1,0), c(1,1), c(0,0))))
p
unclass(p)
a = 1:8
class(a)
attr(a, "dim") = c(2,4) # or: dim(a) = c(2,4)
class(a)
a
attr(a, "dim") = c(2,2,2) # or: dim(a) = c(2,2,2)
class(a)
a
a = c(first = 3, second = 4, last = 5)
a["second"]
attributes(a)
a = matrix(1:4, 2, 2)
dimnames(a) = list(rows = c("row1", "row2"), cols = c("col1", "col2"))
a
attributes(a)
df = data.frame(a = 1:3, b = c(TRUE, FALSE, TRUE))
attributes(df)
f = function(x) {
   a = create_obj(x) # call some other function
   attributes(a) = list(class = "foo", meta = 33)
   a
}
f = function(x) {
   a = create_obj(x) # call some other function
   structure(a, class = "foo", meta = 33)
}
set.seed(1014)
options(digits = 3)

knitr::opts_chunk$set(
  comment = "#",
  collapse = knitr::is_latex_output(),
  cache = TRUE,
#  out.width = "70%",
  fig.align = 'center'
#  fig.width = 6,
#  fig.asp = 0.618,  # 1 / phi
#  fig.show = "hold"
)

options(dplyr.print_min = 6, dplyr.print_max = 6)
options(stars.crs = 17)
mapview::mapviewOptions(fgb = FALSE)

# for units: (?)
Sys.setenv(UDUNITS2_XML_PATH="")
if (knitr::is_latex_output())
    options(width = 66)
    #options(width = 72)