rggplot2r-sfkernel-density

KDE using simple features grid uses centroids instead of grid cell

I would like to make a KDE out of an sf object. For this, I made a grid as rectangle, but R does not take it. Instead, it is using the centroid of the polygons.

This is what I tried:

library(sf)
library(ggplot2)
library(SpatialKDE)
library(rnaturalearth)
library(rnaturalearthdata)

d_trees <- sf::st_read("digitised_trees.shp")
# Output: Reading layer `digitised_trees' from data source 
# `C:\Users\regin\Documents\UNIGIS_Salzburg\Modul_8\aufgabe_2\digitised_trees.shp' using driver `ESRI Shapefile'
# Simple feature collection with 4145 features and 2 fields
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: 672578 ymin: 5188850 xmax: 673748.4 ymax: 5189592
# Projected CRS: WGS 84 / UTM zone 32N

s_trees <- sf::st_read("simulated_trees.shp")
# Output: Reading layer `simulated_trees' from data source 
# `C:\Users\regin\Documents\UNIGIS_Salzburg\Modul_8\aufgabe_2\simulated_trees.shp' using driver `ESRI Shapefile'
# Simple feature collection with 3879 features and 2 fields
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: 672633.5 ymin: 5188850 xmax: 673749.2 ymax: 5189561
# Projected CRS: WGS 84 / UTM zone 32N

d_trees_subset <- d_trees[1:10, ]
s_trees_subset <- s_trees[1:10, ]
dput(d_trees_subset)
structure(list(HEIGHT = c(">15", ">15", ">15", ">15", ">15", 
">15", "6-15", "6-15", "6-15", "6-15"), elev = c(1673L, 1655L, 
1661L, 1716L, 1699L, 1725L, 1741L, 1768L, 1759L, 1762L), geometry = structure(list(
    structure(c(672791.910571324, 5188901.84124404), class = c("XY", 
    "POINT", "sfg")), structure(c(672799.480667902, 5188865.25244391
    ), class = c("XY", "POINT", "sfg")), structure(c(672794.43393685, 
    5188880.39263707), class = c("XY", "POINT", "sfg")), structure(c(672680.146506565, 
    5188881.75946006), class = c("XY", "POINT", "sfg")), structure(c(672694.445577879, 
    5188860.73141401), class = c("XY", "POINT", "sfg")), structure(c(672651.548363936, 
    5188869.14263243), class = c("XY", "POINT", "sfg")), structure(c(672670.315895036, 
    5188914.77349236), class = c("XY", "POINT", "sfg")), structure(c(672662.745798458, 
    5188954.5164994), class = c("XY", "POINT", "sfg")), structure(c(672682.254178773, 
    5188954.91661639), class = c("XY", "POINT", "sfg")), structure(c(672670.946736418, 
    5188948.83892696), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT", 
"sfc"), precision = 0, bbox = structure(c(xmin = 672651.548363936, 
ymin = 5188860.73141401, xmax = 672799.480667902, ymax = 5188954.91661639
), class = "bbox"), crs = structure(list(input = "EPSG:32632", 
    wkt = "PROJCRS[\"WGS 84 / UTM zone 32N\",\n    BASEGEOGCRS[\"WGS 84\",\n        ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n            MEMBER[\"World Geodetic System 1984 (Transit)\"],\n            MEMBER[\"World Geodetic System 1984 (G730)\"],\n            MEMBER[\"World Geodetic System 1984 (G873)\"],\n            MEMBER[\"World Geodetic System 1984 (G1150)\"],\n            MEMBER[\"World Geodetic System 1984 (G1674)\"],\n            MEMBER[\"World Geodetic System 1984 (G1762)\"],\n            MEMBER[\"World Geodetic System 1984 (G2139)\"],\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]],\n            ENSEMBLEACCURACY[2.0]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4326]],\n    CONVERSION[\"UTM zone 32N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",9,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Navigation and medium accuracy spatial referencing.\"],\n        AREA[\"Between 6°E and 12°E, northern hemisphere between equator and 84°N, onshore and offshore. Algeria. Austria. Cameroon. Denmark. Equatorial Guinea. France. Gabon. Germany. Italy. Libya. Liechtenstein. Monaco. Netherlands. Niger. Nigeria. Norway. Sao Tome and Principe. Svalbard. Sweden. Switzerland. Tunisia. Vatican City State.\"],\n        BBOX[0,6,84,12]],\n    ID[\"EPSG\",32632]]"), class = "crs"), n_empty = 0L)), sf_column = "geometry", agr = structure(c(HEIGHT = NA_integer_, 
elev = NA_integer_), levels = c("constant", "aggregate", "identity"
), class = "factor"), row.names = c(NA, 10L), class = c("sf", 
"data.frame"))
dput(s_trees_subset)
structure(list(HEIGHT = c(34.9625800026028, 34.7512513942823, 
34.8683319241762, 34.8568493729999, 34.9311321770029, 34.8433027568252, 
34.9012116552477, 34.9123849337441, 34.9424040445574, 34.9508925408629
), elev = c(1587L, 1607L, 1562L, 1610L, 1652L, 1532L, 1597L, 
1553L, 1528L, 1563L), geometry = structure(list(structure(c(673656.645769361, 
5188907.16103378), class = c("XY", "POINT", "sfg")), structure(c(673025.118431594, 
5188931.72983531), class = c("XY", "POINT", "sfg")), structure(c(673453.162168716, 
5188928.38909054), class = c("XY", "POINT", "sfg")), structure(c(673440.619588101, 
5189028.38824619), class = c("XY", "POINT", "sfg")), structure(c(673589.648795813, 
5189023.50717526), class = c("XY", "POINT", "sfg")), structure(c(673448.33074397, 
5188884.22643243), class = c("XY", "POINT", "sfg")), structure(c(673558.477587886, 
5188944.2819144), class = c("XY", "POINT", "sfg")), structure(c(673140.127848215, 
5188976.35201315), class = c("XY", "POINT", "sfg")), structure(c(673191.602010456, 
5189012.64219437), class = c("XY", "POINT", "sfg")), structure(c(673079.423488968, 
5188931.033843), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT", 
"sfc"), precision = 0, bbox = structure(c(xmin = 673025.118431594, 
ymin = 5188884.22643243, xmax = 673656.645769361, ymax = 5189028.38824619
), class = "bbox"), crs = structure(list(input = "EPSG:32632", 
    wkt = "PROJCRS[\"WGS 84 / UTM zone 32N\",\n    BASEGEOGCRS[\"WGS 84\",\n        ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n            MEMBER[\"World Geodetic System 1984 (Transit)\"],\n            MEMBER[\"World Geodetic System 1984 (G730)\"],\n            MEMBER[\"World Geodetic System 1984 (G873)\"],\n            MEMBER[\"World Geodetic System 1984 (G1150)\"],\n            MEMBER[\"World Geodetic System 1984 (G1674)\"],\n            MEMBER[\"World Geodetic System 1984 (G1762)\"],\n            MEMBER[\"World Geodetic System 1984 (G2139)\"],\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]],\n            ENSEMBLEACCURACY[2.0]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4326]],\n    CONVERSION[\"UTM zone 32N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",9,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Navigation and medium accuracy spatial referencing.\"],\n        AREA[\"Between 6°E and 12°E, northern hemisphere between equator and 84°N, onshore and offshore. Algeria. Austria. Cameroon. Denmark. Equatorial Guinea. France. Gabon. Germany. Italy. Libya. Liechtenstein. Monaco. Netherlands. Niger. Nigeria. Norway. Sao Tome and Principe. Svalbard. Sweden. Switzerland. Tunisia. Vatican City State.\"],\n        BBOX[0,6,84,12]],\n    ID[\"EPSG\",32632]]"), class = "crs"), n_empty = 0L)), sf_column = "geometry", agr = structure(c(HEIGHT = NA_integer_, 
elev = NA_integer_), levels = c("constant", "aggregate", "identity"
), class = "factor"), row.names = c(NA, 10L), class = c("sf", 
"data.frame"))

# Visualising data
ggplot(data = d_trees_subset) +
  geom_sf() + 
  coord_sf(xlim = c(672578, 673748.4), ylim = c(5188850, 5189592)) +
  ggtitle("Digitised Trees")

ggplot(data = s_trees_subset) +
  geom_sf() + 
  coord_sf(xlim = c(672633.5, 673749.2), ylim = c(5188850, 5189561)) +
  ggtitle("Simulated Trees")

band_width <- 50000  
cell_size <- 20000  

# Create grid for KDE
grid_trees <- SpatialKDE::create_grid_rectangular(s_trees_subset, cell_size = cell_size, side_offset = 0)

# Calculate KDE
kde <- SpatialKDE::kde(s_trees_subset, band_width = band_width, kernel = "quartic", grid = grid_trees)
# Output: Using centroids instead of provided `grid` geometries to calculate KDE estimates.

# Load world map
world <- rnaturalearth::ne_countries(scale = 'medium', type = 'map_units', returnclass = 'sf')

# Visualising KDE
ggplot(data = world) +
  geom_sf() +
  geom_sf(data = kde, aes(fill = kde_value), colour = NA) +
  scale_fill_viridis_c(trans = "sqrt", alpha = .6, name = "KDE Value") +
  coord_sf(xlim = c(672578, 673748.4), ylim = c(5188850, 5189592), expand = FALSE) +
  theme(panel.grid.major = element_line(colour = "transparent"))

Can somebody help to resolve this message:

Using centroids instead of provided grid geometries to calculate KDE estimates.

and for the plot, this error:

Error: Invalid index: field name 'x_start' not found
In addition: Warning message:
In st_is_longlat(x) :
bounding box has potentially an invalid value range for longlat data

I did try to change band_width and cell_size a lot, but this did not help at all.

Solution

  • You have to create a raster object in order to produce your desired outcome. To do this, you need to use the raster package. However, raster has been deprecated so your workflow may not work in the future. I suggest looking into alternative KDE methods that are compatible with the terra package to future-proof your analysis.

    As mentioned by @Friede, your error message is from the data in your plot having two different CRS. Your could project your KDE data to match the CRS of your world map, however, there is no point including ne_countries() in your plot. That's because no features will show due to the relatively small extent defined by coord_sf(). In other words, you will not see any distinguishable map features at the extent that you have defined.

    I have included plots using your original method and using a RaterLayer object for comparison.

    library(SpatialKDE)
library(raster)
library(ggplot2)
library(patchwork)

# Set variables
band_width <- 50000  
# Cell size reduced for illustrative purposes due to reduced extent of sampled data
cell_size <- 20

# Your original method
grid_trees <- create_grid_rectangular(s_trees_subset,
                                      cell_size = cell_size,
                                      side_offset = 0)

# Calculate KDE using sf
kde <- kde(s_trees_subset,
           band_width = band_width,
           kernel = "quartic",
           grid = grid_trees)
# Using centroids instead of provided `grid` geometries to calculate KDE estimates.

# Plot KDE using sf
p1 <- ggplot() +
  geom_sf(data = kde, aes(fill = kde_value), colour = NA) +
  scale_fill_viridis_c(trans = "sqrt", alpha = .6, name = "KDE Value") +
  labs(title = "sf") +
  coord_sf(datum = 25832,
           xlim = c(672578, 673748.4),
           ylim = c(5188850, 5189592),
           expand = FALSE) +
  theme(panel.grid.major = element_line(colour = "transparent"))

# Create RasterLayer grid for KDE
grid_trees <- raster(extent(s_trees_subset), 
                      resolution = cell_size, 
                      crs = crs(s_trees_subset)) 

# Calculate KDE using RasterLayer
kde <- kde(s_trees_subset,
           band_width = band_width,
           kernel = "quartic",
           grid = grid_trees)

# Convert KDE RasterLayer to dataframe for plotting
kde <- as.data.frame(kde, xy = TRUE)

# Plot KDE using RasterLayer
p2 <- ggplot() +
  geom_tile(data = kde,
            aes(x, y, fill = layer),
            colour = NA) +
  scale_fill_viridis_c(trans = "sqrt", alpha = .6, name = "KDE Value") +
  labs(title = "RasterLayer") +
  coord_sf(xlim = c(672578, 673748.4),
           ylim = c(5188850, 5189592),
           expand = FALSE) +
  theme(panel.grid.major = element_line(colour = "transparent"),
        axis.title = element_blank())

p1 + p2

    1