Search code examples
rspatial

spatial kernel density estimate but i get this error "Error: [rast] cannot create a raster geometry from a single x coordinate"


library(sf)
library(spatialEco)
library(dplyr)

load("C:/Users/user/Desktop/Movedata/geo_record.RData")

geo_data <- geo_record %>%
  rename(x = coords.x1,
         y = coords.x2)
movement_points <- geo_data[-c(9)]

movement_points$longitude <- movement_points$x
movement_points$latitude <-movement_points$y
movement_points_sf <- st_as_sf(movement_points, coords = c("x", "y"), crs = 4326, 
                  agr = "constant")
pt.kde <- sf.kde(x = movement_points_sf, bw = 5, standardize = TRUE, res=40)

i got an error pt.kde <- sf.kde(x = movement_points_sf, bw = 5, standardize = TRUE, res=40)

calculating unweighted kde

Error: [rast] cannot create a raster geometry from a single x coordinate I am using RStudio


Solution

  • The error you are getting is likely because your data are defined by a geographic coordinate system (4326) and sf.kde() is expecting projected data e.g. 3857. Here is a workflow with some example data to illustrate. This example uses the Pseudo-Mercator (EPSG:3857) projected coordinate system. You will need to ensure that the projection you choose is correct for your data and location.

    Note that because you are using dplyr, you can pipe those separate lines of code:

    library(sf)
    library(spatialEco)
    library(dplyr)
    
    # Create some example data
    set.seed(1)
    geo_record <- data.frame(coords.x1 = runif(100, 45.000, 45.002),
                             coords.x2 = runif(100, 45.000, 45.002))
    
    geo_record <- geo_record %>%
      mutate(!!!setNames(rep(NA, length(letters[1:7])), letters[1:7]))
    
    geo_record <- geo_record %>%
      mutate(across(a:g, ~ replace(.x, is.na(.x), 
                                sample(1:10, sum(is.na(.x)), replace = TRUE))))
    
    head(geo_record)
      coords.x1 coords.x2  a  b c  d e f  g
    1  45.00053  45.00131 10  2 5  7 6 5  6
    2  45.00074  45.00071 10  1 5 10 4 8  5
    3  45.00115  45.00054  2  8 3  5 8 3  5
    4  45.00182  45.00199  2  5 5  8 2 8  7
    5  45.00040  45.00127  3  8 8  5 4 5  9
    6  45.00180  45.00043  1 10 7  3 2 3 10
    
    # Modify column names, remove unwanted column, add new columns, convert to sf
    movement_points_sf <- geo_record %>%
      rename(x = coords.x1,
             y = coords.x2) %>%
      select(-9) %>%
      mutate(longitude = x,
             latitude = y) %>%
      st_as_sf(., 
               coords = c("x", "y"), 
               crs = 4326, 
               agr = "constant") %>%
      st_transform(3857)
    
    head(movement_points_sf)
    # Simple feature collection with 6 features and 8 fields
    # Attribute-geometry relationships: constant (8)
    # Geometry type: POINT
    # Dimension:     XY
    # Bounding box:  xmin: 5009422 ymin: 5621589 xmax: 5009579 ymax: 5621834
    # Projected CRS: WGS 84 / Pseudo-Mercator
       a  b c  d e f longitude latitude                geometry
    1 10  2 5  7 6 5  45.00053 45.00131 POINT (5009436 5621728)
    2 10  1 5 10 4 8  45.00074 45.00071 POINT (5009460 5621633)
    3  2  8 3  5 8 3  45.00115 45.00054 POINT (5009505 5621607)
    4  2  5 5  8 2 8  45.00182 45.00199 POINT (5009579 5621834)
    5  3  8 8  5 4 5  45.00040 45.00127 POINT (5009422 5621721)
    6  1 10 7  3 2 3  45.00180 45.00043 POINT (5009577 5621589)
    
    # Calculate KDE
    pt.kde <- sf.kde(x = movement_points_sf, 
                     bw = 5, 
                     standardize = TRUE, 
                     res = 40)
    
    pt.kde
    class       : SpatRaster 
    dimensions  : 4, 7, 1  (nrow, ncol, nlyr)
    resolution  : 28.57143, 80  (x, y)
    extent      : 5009394, 5009594, 5621486, 5621806  (xmin, xmax, ymin, ymax)
    coord. ref. : WGS 84 / Pseudo-Mercator (EPSG:3857) 
    source(s)   : memory
    name        :        z 
    min value   :    0.000 
    max value   : 3274.323