Create a set of 1000000 (1e6) points, preferably spread (approximately) regularly over the Earth’s surface

O.K., we do a bit less.

library(sf)
## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(stars)
## Loading required package: abind
library(sp)
st_bbox() |>
  st_as_sfc() -> bb
as(bb, "Spatial") |>
  spsample(n = 1000, type = "Fibonacci") |> # https://arxiv.org/pdf/0912.4540.pdf
  st_as_sfc() -> x
## Warning in proj4string(obj): CRS object has comment, which is lost in output; in tests, see
## https://cran.r-project.org/web/packages/sp/vignettes/CRS_warnings.html
plot(x)
plot(bb, add = TRUE)

Or, in orthographic:

par(mfrow=c(2,2), mar = rep(0, 4))
plot(st_transform(x, st_crs("+proj=ortho +lon_0=45 +lat_0=-45")), pch = 3)
plot(st_transform(x, st_crs("+proj=ortho +lon_0=0 +lat_0=-45")), pch = 3)
plot(st_transform(x, st_crs("+proj=ortho +lon_0=-30 +lat_0=35")), pch = 3)
plot(st_transform(x, st_crs("+proj=ortho +lon_0=-90 +lat_0=-15")), pch = 3)

Some more points:

as(bb, "Spatial") |>
  spsample(n = 1e5, type = "Fibonacci") |>
  st_as_sfc() -> x
## Warning in proj4string(obj): CRS object has comment, which is lost in output; in tests, see
## https://cran.r-project.org/web/packages/sp/vignettes/CRS_warnings.html

Subset this dataset for points that are on land

library(rnaturalearth)
co = ne_countries(returnclass = "sf") |> st_make_valid()
x = x[co]
plot(x)

Generate for the entire Earth, …

… on a regular grid, densities of these land points (also called intensity: number of points per unit area), using different bandwidths (e.g. 50, 100, 200, 500 and 1000 km)

Create a function that takes a point set and a bandwidth as input arguments and that returns a global grid with point densities

(grd = st_as_stars(st_bbox(), nx = 90, ny = 45)) # creates 4-degree global grid
## stars object with 2 dimensions and 1 attribute
## attribute(s):
##         Min. 1st Qu. Median Mean 3rd Qu. Max.
## values     0       0      0    0       0    0
## dimension(s):
##   from to offset delta refsys point values x/y
## x    1 90   -180     4 WGS 84    NA   NULL [x]
## y    1 45     90    -4 WGS 84    NA   NULL [y]

Slow:

calc_dens = function(pts, grd, dist) {
        buf = st_buffer(pts, dist)
        a = st_area(buf) |> units::set_units(km^2) |> units::drop_units()
        lengths(st_intersects(grd, buf)) / a
}
grd$d1000 = calc_dens(x, st_as_sfc(grd, as_points = TRUE), units::set_units(1000, km))
plot(grd["d1000"])
title("points per square km")

Faster:

calc_dens = function(pts, grd, dist) {
        pt = st_sfc(st_point(c(0,0)), crs = 4326)
        a = st_area(st_buffer(pt, dist)) |> units::set_units(km^2) |> units::drop_units()
        lengths(st_is_within_distance(grd, pts, dist)) / a
}
grd$d100 = calc_dens(x, st_as_sfc(grd, as_points = TRUE), units::set_units(100, km))
grd$d200 = calc_dens(x, st_as_sfc(grd, as_points = TRUE), units::set_units(200, km))
grd$d500 = calc_dens(x, st_as_sfc(grd, as_points = TRUE), units::set_units(500, km))
grd$d1000 = calc_dens(x, st_as_sfc(grd, as_points = TRUE), units::set_units(1000, km))
grd$values = NULL
plot(merge(grd))

Visualise these maps with densities using an orthographic projection, from several different perspectives

par(mfrow=c(2,2), mar = c(0, 0, 1.1, 0))
sf = st_as_sf(grd)
ne = function(x, lat, lon) { 
        crs = st_crs(paste0("+proj=ortho +lon_0=", lon, " +lat_0=", lat))
        pt = st_sfc(st_point(c(lon, lat)), crs = 4326)
        bu = st_buffer(pt, units::set_units(9800000, m))
        x = st_transform(st_intersection(x, bu), crs);
        plot(x, key.pos = NULL, border = NA, reset = FALSE)
        plot(st_transform(st_intersection(st_geometry(co), bu), crs), add = TRUE)
}
ne(sf[4], 45, -45)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
ne(sf[4], -30, 45)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
ne(sf[4], -45, -135)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
ne(sf[4], 90, 0)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

Create a spinning globe (animated) plot of the densities

(for bonus points) build in the option to choose different kernel functions for density computation

(for bonus points) let the densities not spill over in water areas, i.e. compute them only over land and take edge effects into account

This only clips the maps above, it doesn’t compute taking water into account:

plot(merge(grd)[co])

Spherical Geometry