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",
  "dbscan",
  "dbscan",
  "ecmwfr",
  "gdalcubes",
  "gdalUtils",
  "ggplot2",
  "ggspatial",
  "gstat",
  "hglm",
  "HSAR",
  "igraph",
  "INLA",
  "lme4",
  "lwgeom",
  "mapsf",
  "mapview",
  "Matrix",
  "mgcv",
  "osmar", 
  "R2BayesX",
  "raster",
  "RColorBrewer",
  "rgee",
  "rgeoda",
  "rnaturalearth",
  "rstac",
  "s2",
  "sfnetworks",
  "sf", 
  "sp",
  "spacetime",
  "spatialreg",
  "spatstat",
  "spData",
  "spdep",
  "splm",
  "stars",
  "stcos",
  "stplanr", 
  "terra",
  "tidyverse",
  "tmap",
  "tsibble",
  "units",
  "viridis"

  ), "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, 'EPSG: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) %>% 
    pivot_longer(starts_with("SID")) -> nc_longer
ggplot() + geom_sf(data = nc_longer, aes(fill = value)) + 
  facet_wrap(~ name, ncol = 1, labeller = labeller(name = 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))
mapviewOptions(fgb = FALSE)
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')
s2_orig = 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:
sf_use_s2(s2_orig) # back to what it was before changing
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, 3))
set.seed(133331)
mp = st_multipoint(matrix(runif(20), 10))
plot(mp, cex = 2)
plot(st_convex_hull(mp), add = TRUE, col = NA, border = 'red')
box()
plot(mp, cex = 2)
plot(st_voronoi(mp), add = TRUE, col = NA, border = 'red')
box()
plot(mp, cex = 2)
plot(st_triangulate(mp), add = TRUE, col = NA, border = 'darkgreen')
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 = 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, breaks = "equal", main = NA) # 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"))
# 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))
g = st_make_grid(st_bbox(st_as_sfc("LINESTRING(0 0,1 1)")), n = c(2,2))
par(mar = rep(0,4))
plot(g)
plot(g[1] * diag(c(3/4, 1)) + c(0.25, 0.125), add = TRUE, lty = 2)
text(c(.2, .8, .2, .8), c(.2, .2, .8, .8), c(1,2,4,8), col = 'red')
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/master/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(st_geometry(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)
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))
set.seed(131)
circles = st_sample(st_as_sfc(st_bbox(r)), 3) %>%
    st_buffer(500)
aggregate(r, circles, max)
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_drop_geometry(trn))
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) %>% 
    pivot_longer(3:6, names_to = "mode", values_to = "n")
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(adrop(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)
(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),]
cat(st_crs("EPSG:4326")$wkt, "\n")
st_axis_order(TRUE)
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) &&
    !inherits(try(DBI::dbConnect( 
        RPostgres::Postgres(),
        host = "localhost",
        dbname = "postgis"), silent = TRUE), "try-error")
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)
options(timeout = 600) # or large in case of slow network
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)
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
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) %>% 
    pivot_longer(starts_with("SID")) -> nc_longer
ggplot() + geom_sf(data = nc_longer, aes(fill = value)) + 
  facet_wrap(~ name, ncol = 1, labeller = labeller(name = year_labels)) +
  scale_y_continuous(breaks = 34:36) +
  scale_fill_gradientn(colours = sf.colors(20)) +
  theme(panel.grid.major = element_line(colour = "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)
system.file("gpkg/nc.gpkg", package = "sf") %>%
    read_sf() %>%
    st_transform('EPSG:32119') -> nc.32119
tm_shape(nc.32119) + tm_polygons(c("SID74", "SID79"))
nc.32119 %>% select(SID74, SID79) %>% 
    pivot_longer(starts_with("SID"), values_to = "SID") -> nc_longer
tm_shape(nc_longer) + tm_polygons("SID") + tm_facets(by = "name")
nc.32119 %>% select(SID74, SID79) %>% 
    pivot_longer(starts_with("SID"), values_to = "SID") -> nc_longer
tm_shape(nc_longer) + tm_polygons("SID") + tm_facets(by = "name")
tm_shape(r) + tm_raster()
tm_shape(r) + tm_raster()
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(13531)
library(sf)
n = 30
xy = data.frame(x = runif(n), y = runif(n)) %>% st_as_sf(coords = c("x", "y"))
w1 = st_bbox(c(xmin = 0, ymin = 0, xmax = 1, ymax = 1)) %>%
        st_as_sfc() 
w2 = st_sfc(st_point(c(1, 0.5))) %>% st_buffer(1.2)
par(mfrow = c(1, 2), mar = c(2.1, 2.1, 0.1, 0.5), xaxs = "i", yaxs = "i")
plot(w1, axes = TRUE, col = 'grey')
plot(xy, add = TRUE)
plot(w2, axes = TRUE, col = 'grey')
plot(xy, add = TRUE, cex = .5)
suppressPackageStartupMessages(library(spatstat))
as.ppp(xy)
(pp1 = c(w1, st_geometry(xy)) %>% as.ppp())
c1 = st_buffer(st_centroid(w2), 1.2)
(pp2 = c(c1, st_geometry(xy)) %>% as.ppp())
par(mfrow = c(1, 2), mar = rep(0, 4))
q1 = quadratcount(pp1, nx=3, ny=3)
q2 = quadratcount(pp2, nx=3, ny=3)
plot(q1, main = "")
plot(xy, add = TRUE)
plot(q2, main = "")
plot(xy, add = TRUE)
quadrat.test(pp1, nx=3, ny=3)
quadrat.test(pp2, nx=3, ny=3)
den1 <- density(pp1, sigma = bw.diggle)
den2 <- density(pp2, sigma = bw.diggle)
par(mfrow = c(1, 2), mar = c(0,0,1.1,2))
plot(den1)
plot(pp1, add=TRUE)
plot(den2)
plot(pp1, add=TRUE)
library(stars)
s1 = st_as_stars(den1)
(s2 = st_as_stars(den2))
sum(s1[[1]], na.rm = TRUE)*st_dimensions(s1)$x$delta^2
sum(s2[[1]], na.rm = TRUE)*st_dimensions(s2)$x$delta^2
pt = st_sfc(st_point(c(0.5, 0.5)))
s2$dist = st_distance(st_as_sf(s2, as_points = TRUE, na.rm = FALSE), pt)
(m = ppm(pp2 ~ dist, data = list(dist = as.im(s2["dist"]))))
plot(m, se = FALSE)
predict(m, covariates = list(dist = as.im(s2["dist"]))) %>%
    st_as_stars()
system.file("gpkg/nc.gpkg", package = "sf") %>% 
    read_sf() %>%
    st_geometry() %>%
    st_centroid() %>%
    as.ppp()
longleaf
ll = st_as_sf(longleaf)
print(ll, n = 5)
as.ppp(ll)
print(st_as_sf(copper$SouthLines), n = 5)
print(st_as_sf(chicago), n = 5)
table(st_as_sf(chicago)$label)
kappa = 30 / st_area(w2) # intensity
th = st_sample(w2, kappa = kappa, mu = 3, scale = 0.05, type = "Thomas")
nrow(th)
par(mar = rep(0, 4))
plot(w2)
plot(th, add = TRUE)
# example from plotting chapter:
par(mar = rep(0, 4))
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, border = 'grey')
# sampling from globe:
#sf_use_s2(FALSE)
assign(".sf.use_s2", FALSE, envir=sf:::.sf_cache) # cheat!
pts = suppressMessages( # cheat!
   st_sample(st_as_sfc(st_bbox(st_as_stars())), 1000, exact = FALSE))
#sf_use_s2(TRUE)
assign(".sf.use_s2", TRUE, envir = sf:::.sf_cache) # cheat!
pts = s2_intersection(i, pts) %>% st_as_sfc()
plot(st_transform(pts, "+proj=ortho +lat_0=-10 +lon_0=-30"), 
     add = TRUE, pch = 3, cex = .7)
# we need to detach, to avoid name clashes of using idw() and diagnose() later on
library(spatstat)
w = function(x) which(x == search())
detach(pos = w("package:spatstat"))
detach(pos = w("package:spatstat.linnet"))
detach(pos = w("package:spatstat.core"))
detach(pos = w("package:spatstat.geom"))
detach(pos = w("package:spatstat.data"))
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)
# after running the "Spatiotemporal Geostatistics" chapter, write NO2 means by
write_csv(right_join(a2, tb), "no2.csv")
library(tidyverse)
no2 = read_csv(system.file("external/no2.csv", package = "gstat"))
library(sf)
crs = st_crs("EPSG:32632")
no2.sf = st_as_sf(no2, coords = c("station_longitude_deg", "station_latitude_deg"), crs = "OGC:CRS84") %>%
    st_transform(crs)
data(air, package = "spacetime") # this loads German boundaries into DE_NUTS1
de <- st_transform(st_as_sf(DE_NUTS1), crs)
ggplot() + geom_sf(data = de) +  geom_sf(data = no2.sf, mapping = aes(col = NO2))
# build a grid over Germany:
library(stars)
st_bbox(de) %>%
  st_as_stars(dx = 10000) %>%
  st_crop(de) -> grd
grd
library(gstat)
i = idw(NO2~1, no2.sf, grd)
ggplot() + geom_stars(data = i, aes(fill = var1.pred, x = x, y = y)) + 
    geom_sf(data = st_cast(de, "MULTILINESTRING")) + 
    geom_sf(data = no2.sf)
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)
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) +
    coord_sf(lims_method = "geometry_bbox")
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) %>% 
        pivot_longer(1:2, names_to = "var", values_to = "NO2") -> b2
b2$var = factor(b2$var, levels = c("sample", "kriging"))
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) %>% 
        pivot_longer(1:2, names_to = "var", values_to = "Standard_error") -> b2
b2$var = factor(b2$var, levels = c("sample", "kriging"))
ggplot() + geom_sf(data = b2, mapping = aes(fill = Standard_error)) + facet_wrap(~var, as.table = FALSE) +
     scale_fill_gradientn(colors = sf.colors(20))
s = krige(NO2~1, no2.sf, grd, v.m, nmax = 30, nsim = 6)
library(viridis)
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) +
    coord_sf(lims_method = "geometry_bbox")
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"))
set.seed(123)
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 in local time zone
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)
sel = apply(aq, 2, function(x) mean(is.na(x)) < 0.25)
aqsel = aq[, sel]
library(tidyverse)
read.csv("aq/AirBase_v8_stations.csv", sep = "\t") %>% 
    as_tibble()  %>% 
    filter(country_iso_code == "DE", station_type_of_area == "rural", 
                 type_of_station == "Background") -> a2
library(sf)
a2.sf = st_as_sf(a2, coords = c("station_longitude_deg", "station_latitude_deg"), crs = 'EPSG: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 = "EPSG: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)
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)]
library(stars)
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))
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(~as.Date(time)) +
    geom_sf(data = st_cast(de, "MULTILINESTRING")) + 
    geom_sf(data = no2.sf, col = 'grey', cex = .5) + 
    coord_sf(lims_method = "geometry_bbox")
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::opts_chunk$set(echo = TRUE, paged.print=FALSE)
library(sf)
data(pol_pres15, package = "spDataLarge")
pol_pres15 %>% 
    subset(select = c(TERYT, name, types)) %>% 
    head()
library(tmap)
tm_shape(pol_pres15) + tm_fill("types")
if (!all(st_is_valid(pol_pres15))) pol_pres15 <- st_make_valid(pol_pres15)
library(spdep)
args(poly2nb)
system.time(pol_pres15 %>% poly2nb(queen = TRUE) -> nb_q)
nb_q
old_use_s2 <- sf_use_s2()
sf_use_s2(TRUE)
(pol_pres15 %>% st_transform("OGC:CRS84") -> pol_pres15_ll) %>% 
    poly2nb(queen = TRUE) -> nb_q_s2
all.equal(nb_q, nb_q_s2, check.attributes=FALSE)
(nb_q %>% n.comp.nb())$nc
library(Matrix, warn.conflicts = FALSE)
nb_q %>% 
    nb2listw(style = "B") %>% 
    as("CsparseMatrix") -> smat
library(igraph, warn.conflicts = FALSE)
(smat %>% 
        graph.adjacency() -> g1) %>% 
    count_components()
tf <- tempfile(fileext = ".gal")
write.nb.gal(nb_q, tf)
pol_pres15 %>% 
    st_geometry() %>% 
    st_centroid(of_largest_polygon = TRUE) -> coords 
(coords %>% tri2nb() -> nb_tri)
nb_tri %>% 
    nbdists(coords) %>% 
    unlist() %>% 
    summary()
(nb_tri %>% n.comp.nb())$nc
(nb_tri %>% 
        soi.graph(coords) %>% 
        graph2nb() -> nb_soi)
(nb_soi %>% n.comp.nb() -> n_comp)$nc
table(n_comp$comp.id)
opar <- par(mar = c(0,0,0,0)+0.5)
pol_pres15 %>% 
    st_geometry() %>% 
    plot(border = "grey", lwd = 0.5)
nb_soi %>% 
    plot(coords = coords, add = TRUE, points = FALSE, lwd = 0.5)
nb_tri %>% 
    diffnb(nb_soi) %>% 
    plot(coords = coords, col = "orange", add = TRUE, points = FALSE, lwd = 0.5)
par(opar)
coords %>% 
    knearneigh(k = 1) %>% 
    knn2nb() %>% 
    nbdists(coords) %>% 
    unlist() %>% 
    summary()
system.time(coords %>% dnearneigh(0, 18000) -> nb_d18)
system.time(coords %>% dnearneigh(0, 18000, use_kd_tree = FALSE) -> nb_d18a)
all.equal(nb_d18, nb_d18a, check.attributes = FALSE)
nb_d18
(nb_d18 %>% n.comp.nb() -> n_comp)$nc
table(n_comp$comp.id)
(coords %>% dnearneigh(0, 18300) -> nb_d183)
(nb_d183 %>% n.comp.nb())$nc
(coords %>% dnearneigh(0, 16000) -> nb_d16)
((coords %>% knearneigh(k = 6) -> knn_k6) %>% knn2nb() -> nb_k6)
(knn_k6 %>% knn2nb(sym = TRUE) -> nb_k6s)
(nb_k6s %>% n.comp.nb())$nc
old_use_s2 <- sf_use_s2()
sf_use_s2(TRUE)
pol_pres15_ll %>% 
    st_geometry() %>% 
    st_centroid(of_largest_polygon = TRUE) -> coords_ll
(coords_ll %>% knearneigh(k = 6) %>% knn2nb() -> nb_k6_ll)
isTRUE(all.equal(nb_k6, nb_k6_ll, check.attributes = FALSE))
nb_q %>% nbdists(coords_ll) %>% unlist() %>% summary()
nb_q %>% nbdists(coords) %>% unlist() %>% summary()
sf_use_s2(old_use_s2)
args(nb2listw)
args(spweights.constants)
(nb_q %>% 
        nb2listw(style = "B") -> lw_q_B) %>% 
    spweights.constants() %>% 
    data.frame() %>% 
    subset(select = c(n, S0, S1, S2))
(nb_q %>% 
        nb2listw(style = "W") -> lw_q_W) %>% 
    spweights.constants() %>% 
    data.frame() %>% 
    subset(select = c(n, S0, S1, S2))
nb_d183 %>% 
    nbdists(coords) %>% 
    lapply(function(x) 1/(x/1000)) -> gwts
(nb_d183 %>% nb2listw(glist=gwts, style="B") -> lw_d183_idw_B) %>% 
    spweights.constants() %>% 
    data.frame() %>% 
    subset(select=c(n, S0, S1, S2))
try(nb_d16 %>% nb2listw(style="B") -> lw_d16_B)
nb_d16 %>% 
    nb2listw(style="B", zero.policy=TRUE) %>% 
    spweights.constants(zero.policy=TRUE) %>% 
    data.frame() %>% 
    subset(select=c(n, S0, S1, S2))
nb_q
(nb_q %>% nblag(2) -> nb_q2)[[2]]
nblag_cumul(nb_q2)
union.nb(nb_q2[[2]], nb_q2[[1]])
diameter(g1)
g1 %>% shortest.paths() -> sps
(sps %>% apply(2, max) -> spmax) %>% max()
mr <- which.max(spmax)
pol_pres15$name0[mr]
pol_pres15$sps1 <- sps[,mr]
tm1 <- tm_shape(pol_pres15) + tm_fill("sps1", title = "Shortest path\ncount")
coords[mr] %>% 
    st_distance(coords) %>% 
    c() %>% 
    (function(x) x/1000)() %>% 
    units::set_units(NULL) -> pol_pres15$dist_52
library(ggplot2)
g1 <- ggplot(pol_pres15, aes(x = sps1, y = dist_52)) + geom_point() + xlab("Shortest path count") +  ylab("km distance")
gridExtra::grid.arrange(tmap_grob(tm1), g1, nrow=1)
glance_htest <- function(ht) c(ht$estimate, 
    "Std deviate" = unname(ht$statistic), 
    "p.value" = unname(ht$p.value))
set.seed(1)
(pol_pres15 %>% 
        nrow() %>% 
        rnorm() -> x) %>% 
    moran.test(lw_q_B, randomisation = FALSE, alternative = "two.sided") %>% 
    glance_htest()
beta <- 0.0015
coords %>% 
    st_coordinates() %>% 
    subset(select = 1, drop = TRUE) %>% 
    (function(x) x/1000)() -> t
(x + beta * t -> x_t) %>% 
    moran.test(lw_q_B, randomisation = FALSE, alternative = "two.sided") %>% 
    glance_htest()
lm(x_t ~ t) %>% 
    lm.morantest(lw_q_B, alternative = "two.sided") %>% 
    glance_htest()
args(joincount.test)
(pol_pres15 %>% 
        st_drop_geometry() %>% 
        subset(select = types, drop = TRUE) -> Types) %>% 
    table()
Types %>% joincount.multi(listw = lw_q_B)
Types %>% joincount.multi(listw = lw_d183_idw_B)
args(moran.test)
(pol_pres15 %>% 
        st_drop_geometry() %>% 
        subset(select = I_turnout, drop = TRUE) -> z) %>% 
    moran.test(listw = lw_q_B, randomisation = FALSE) %>% 
    glance_htest()
lm(I_turnout ~ 1, pol_pres15) %>% 
    lm.morantest(listw = lw_q_B) %>% 
    glance_htest()
(z %>% 
    moran.test(listw = lw_q_B) -> mtr) %>% 
    glance_htest()
set.seed(1)
z %>% 
    moran.mc(listw = lw_q_B, nsim = 999, return_boot = TRUE) -> mmc
c("Permutation bootstrap" = var(mmc$t), 
  "Analytical randomisation" = unname(mtr$estimate[3]))
z %>% 
    moran.plot(listw = lw_q_W, labels = pol_pres15$TERYT, cex = 1, pch = ".",
        xlab = "I round turnout", ylab = "lagged turnout") -> infl_W
pol_pres15$hat_value <- infl_W$hat
tm_shape(pol_pres15) + tm_fill("hat_value")
args(localmoran)
z %>% 
    localmoran(listw = lw_q_W, alternative = "two.sided") -> locm
all.equal(sum(locm[,1])/Szero(lw_q_W), unname(moran.test(z, lw_q_W)$estimate[1]))
pva <- function(pv) cbind("none" = pv, "bonferroni" = p.adjust(pv, "bonferroni"),
    "fdr" = p.adjust(pv, "fdr"), "BY" = p.adjust(pv, "BY"))
locm %>% 
    subset(select = "Pr(z != 0)", drop = TRUE) %>% 
    pva() -> pvsp
f <- function(x) sum(x < 0.05)
apply(pvsp, 2, f)
library(parallel)
set.coresOption(ifelse(detectCores() == 1, 1, detectCores()-1L))
system.time(z %>% 
        localmoran_perm(listw = lw_q_W, nsim = 499, alternative = "two.sided",
            iseed = 1) -> locm_p)
locm_p %>% 
    subset(select = "Pr(z != 0)", drop = TRUE) %>% 
    pva() -> pvsp
apply(pvsp, 2, f)
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))
(locm_p %>% 
        subset(select = Z.Ii, drop = TRUE) %>% 
        cut(brks) %>% 
        table()-> tab)
sum(tab[c(1:5, 8:12)])
sum(tab[c(1, 12)])
z %>% 
    localmoran(listw = lw_q_W, conditional = TRUE, alternative = "two.sided") -> locm_c
locm_c %>% 
    subset(select = "Pr(z != 0)", drop = TRUE) %>% 
    pva() -> pvsp
apply(pvsp, 2, f)
pol_pres15$locm_Z <- locm[, "Z.Ii"]
pol_pres15$locm_c_Z <- locm_c[, "Z.Ii"]
pol_pres15$locm_p_Z <- locm_p[, "Z.Ii"]
tm_shape(pol_pres15) + tm_fill(c("locm_Z", "locm_c_Z", "locm_p_Z"), breaks = brks,
    midpoint = 0, title = "Standard deviates of\nLocal Moran's I") +
    tm_facets(free.scales = FALSE, ncol = 2) + tm_layout(panel.labels = c("Analytical total",
        "Analytical conditional", "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_ac_q <- quadr
is.na(pol_pres15$hs_ac_q) <- !(pol_pres15$locm_c_Z < brks[2] | 
        pol_pres15$locm_c_Z > brks[12])
c <- table(pol_pres15$hs_ac_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])
d <- table(pol_pres15$hs_cp_q)
t(rbind("Moran plot quadrants" = a, "Unadjusted analytical total" = b, 
    "Bonferroni analytical cond." = c, "Bonferroni cond. perm." = d))
tm_shape(pol_pres15) + tm_fill(c("hs_an_q", "hs_ac_q", "hs_cp_q"), colorNA="grey95",
    textNA="Not significant", title="Turnout hotspot status\nLocal Moran's I") +
    tm_facets(free.scales=FALSE, ncol=2) + tm_layout(panel.labels=
            c("Unadjusted analytical total", "Bonferroni analytical cond.", 
              "Cond. perm. with Bonferroni"))
lm(z ~ 1, weights = pol_pres15$I_entitled_to_vote) -> lm_null
system.time(lm_null %>% localmoran.sad(nb = nb_q, style = "W", alternative = "two.sided") %>%
        summary() -> locm_sad_null)
lm(z ~ Types, weights=pol_pres15$I_entitled_to_vote) -> lm_types
system.time(lm_types %>% localmoran.sad(nb = nb_q, style = "W", alternative = "two.sided") %>%
        summary() -> locm_sad_types)
spatialreg::lmSLX(z ~ Types, listw = lw_q_W, 
    weights = pol_pres15$I_entitled_to_vote) -> lm_Dtypes
system.time(lm_Dtypes %>% localmoran.sad(nb = nb_q, style = "W", alternative = "two.sided") %>%
        summary() -> locm_sad_Dtypes)
pol_pres15$locm_sad_null_Z <- locm_sad_null[, "Saddlepoint"]
pol_pres15$locm_sad_types_Z <- locm_sad_types[, "Saddlepoint"]
pol_pres15$locm_sad_Dtypes_Z <- locm_sad_Dtypes[, "Saddlepoint"]
tm_shape(pol_pres15) + tm_fill(c("locm_sad_null_Z", "locm_sad_types_Z",
    "locm_sad_Dtypes_Z", "locm_c_Z"), breaks = brks, midpoint=0,
    title="Standard deviates of\nLocal Moran's I") +
    tm_facets(free.scales = FALSE, ncol = 2) + tm_layout(panel.labels = c(
        "Saddlepoint weighted null",
        "Saddlepoint weighted types",
        "Saddlepoint weighted Durbin types",
        "Analytical conditional"))
system.time(z %>% 
        localG(lw_q_W) -> locG)
system.time(z %>% 
        localG_perm(lw_q_W, nsim = 499, iseed = 1) -> locG_p)
locG %>% 
    c() %>% 
    abs() %>% 
    pnorm(lower.tail = FALSE) %>% 
    (function(x) x*2)() %>% 
    pva() -> pvsp
apply(pvsp, 2, f)
library(ggplot2)
p1 <- ggplot(data.frame(Zi=locm_c[,4], Zi_perm=locm_p[,4])) + geom_point(aes(x=Zi,
    y=Zi_perm), alpha=0.2) + xlab("Analytical conditional") + 
    ylab("Bootstrap conditional") + 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 conditional") + 
    ylab("Bootstrap conditional") + 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="Local G Bonferroni\nadjusted hotspot 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)
knitr::opts_chunk$set(echo = TRUE, paged.print = FALSE)
library(sf)
library(spData)
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_506$CHAS <- as.factor(boston_506$CHAS)
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")
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)))
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)))
library(Matrix)
library(lme4)
MLM <- lmer(update(form, . ~ . + (1 | NOX_ID)), data = boston_487, REML = FALSE)
boston_93$MLM_re <- ranef(MLM)[[1]][,1]
suppressPackageStartupMessages(library(hglm))
suppressWarnings(HGLM_iid <- hglm(fixed = form, random = ~1 | NOX_ID, data = boston_487,
    family = gaussian()))
boston_93$HGLM_re <- unname(HGLM_iid$ranef)
W <- as(spdep::nb2listw(nb_q_93, style = "B"), "CsparseMatrix")
suppressWarnings(HGLM_car <- hglm(fixed = form, random = ~ 1 | NOX_ID, data = boston_487,
    family = gaussian(), rand.family = CAR(D=W)))
boston_93$HGLM_ss <- HGLM_car$ranef[,1]
library(HSAR)
library(MatrixModels, warn.conflicts = FALSE)
Z <- as(model.Matrix(~ -1 + as.factor(NOX_ID), data = boston_487, sparse = TRUE),
    "dgCMatrix")
set.seed(123)
suppressWarnings(HSAR_ss <- hsar(form, data=boston_487, W = NULL, M = W, Delta = Z, 
                              burnin = 2000, Nsim = 12000, thinning = 2))
boston_93$HSAR_ss <- HSAR_ss$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)
all.equal(row.names(RBX_gra), attr(nb_q_93, "region.id"))
all.equal(unname(diag(RBX_gra)), spdep::card(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
suppressPackageStartupMessages(library(INLA))
INLA_iid <- inla(update(form, . ~ . + f(NOX_ID, model = "iid")), family = "gaussian", 
    data = boston_487)
boston_93$INLA_re <- INLA_iid$summary.random$NOX_ID$mean
ID2 <- as.integer(as.factor(boston_487$NOX_ID))
INLA_ss <- inla(update(form, . ~ . + f(ID2, model = "besag", graph = W)), family = "gaussian",
    data = boston_487)
boston_93$INLA_ss <- INLA_ss$summary.random$ID2$mean
M <- Diagonal(nrow(W), rowSums(W)) - W
Cmatrix <- Diagonal(nrow(M), 1) -  M
INLA_lr <- inla(update(form, . ~ . + f(ID2, model = "generic1", Cmatrix = Cmatrix)),
    family = "gaussian", data = boston_487)
boston_93$INLA_lr <- INLA_lr$summary.random$ID2$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")
ssre <- predict(GAM_MRF, type = "terms", se = FALSE)[, "s(NOX_ID)"]
all(sapply(tapply(ssre, list(boston_487$NOX_ID), c), function(x) length(unique(x)) == 1))
boston_93$GAM_ss <- aggregate(ssre, list(boston_487$NOX_ID), head, n=1)$x
library(tmap, warn.conflicts=FALSE)
tm_shape(boston_93) + tm_fill(c("MLM_re", "HGLM_re", "INLA_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("lmer", "hglm", "inla", "bayesx"))
tm_shape(boston_93) + tm_fill(c("HGLM_ss", "HSAR_ss", "INLA_lr", "INLA_ss", "BX_ss",
    "GAM_ss"), midpoint = 0, title = "SSRE")  + tm_facets(free.scales = FALSE) + 
    tm_borders(lwd = 0.3, alpha = 0.4) + tm_layout(panel.labels = c("hglm CAR", "hsar SAR",
        "inla Leroux", "inla ICAR", "bayesx ICAR", "gam ICAR"))
library(spatialreg)
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]
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))
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_489w <- lmSLX(form, data = boston_489, listw = lw_q_489, weights = units, zero.policy = TRUE)
SDEM_489w <- errorsarlm(form, data = boston_489, listw = lw_q_489, Durbin  =TRUE,
    weights = units, zero.policy = TRUE, control = list(pre_eig = eigs_489))
broom::tidy(lmtest::lrtest(SLX_489w, SDEM_489w))
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))
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_SLX <- summary(impacts(SLX_94))
rbind(Impacts = sum_imp_94_SLX$mat[5,], SE = sum_imp_94_SLX$semat[5,])
sum_imp_94_SDEMw <- summary(impacts(SDEM_94w))
rbind(Impacts = sum_imp_94_SDEMw$mat[5,], SE = sum_imp_94_SDEMw$semat[5,])
sum_imp_94_SLXw <- summary(impacts(SLX_94w))
rbind(Impacts = sum_imp_94_SLXw$mat[5,], SE = sum_imp_94_SLXw$semat[5,])
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)))
library(INLA)
W <- as(lw_q, "CsparseMatrix")
n <- nrow(W)
e <- eigenw(lw_q)
re.idx <- which(abs(Im(e)) < 1e-6)
rho.max <- 1/max(Re(e[re.idx]))
rho.min <- 1/min(Re(e[re.idx]))
rho <- mean(c(rho.min, rho.max))
boston_506$idx <- 1:n
zero.variance = list(prec=list(initial = 25, fixed = TRUE))
args.slm <- list(rho.min = rho.min, rho.max = rho.max, W = W, X = matrix(0, n, 0),
    Q.beta = matrix(1,0,0))
hyper.slm <- list(prec = list(prior = "loggamma", param = c(0.01, 0.01)),
    rho = list(initial = 0, prior = "logitbeta", param = c(1,1)))
WX <- create_WX(model.matrix(update(form, CMEDV ~ .), data=boston_506), lw_q)
SDEM_506_slm <- inla(update(form, . ~ . + WX + f(idx, model = "slm", args.slm = args.slm,
    hyper = hyper.slm)), data = boston_506, family = "gaussian", control.family = 
        list(hyper = zero.variance), control.compute = list(dic = TRUE, cpo = TRUE))
mvs <- SDEM_506_slm$marginals.fitted.values[which(is.na(boston_506$median))]
mv_mean <- sapply(mvs, function(x) mean(exp(x[, 1]), ))
mv_sd <- sapply(mvs, function(x) sd(exp(x[, 1])))
data.frame(fit_TS = t0[,1], fit_KP2 = c(t1), fit_KP5 = c(t2), INLA_slm = mv_mean,
    INLA_slm_sd = mv_sd, censored = boston_506$censored[as.integer(attr(t0, "region.id"))])
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(sp)
y = as(x, "Spatial")
x0 = st_as_sf(y)
x = as(sf::read_sf("file"), "Spatial")
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)
}
library(sf)
system.file("gpkg/nc.gpkg", package = "sf") %>%
    read_sf() -> nc
attributes(nc)
nc$geom
attributes(nc$geom)
nc$geom[[4]]
typeof(nc$geom[[4]])
attributes(nc$geom[[4]])
length(nc$geom[[4]])
lengths(nc$geom)
typeof(nc$geom[[4]])
typeof(nc$geom[[4]][[1]])
length(nc$geom[[4]][[1]])
typeof(nc$geom[[4]][[1]][[1]])
dim(nc$geom[[4]][[1]][[1]])
head(nc$geom[[4]][[1]][[1]])
nc$geom[[4]][[1]][[1]][3,2] = 36.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)