Search code examples
rr-sfspatstat

Convert sf to unmarked ppp


I wish to convert an sf object to an unmarked ppp. (Conversion from sf to ppp is now supported, according to this post.)

library(sf)

#Initialise sf object
pp <- structure(list(X = c(959207.877070254, 959660.734838225, 951483.685462513,  951527.767554883, 958310.673042469, 950492.05212104, 959660.734838225,  959207.877070254, 960500.020456073, 959660.734838225), Y = c(1944457.42827898,  1955543.76027363, 1939982.16629396, 1940216.55143212, 1954704.68186897,  1951434.68524296, 1955543.76027363, 1944457.42827898, 1955292.64874361,  1955543.76027363), geometry = structure(list(structure(c(959207.877070254,  1944457.42827898), class = c("XY", "POINT", "sfg")), structure(c(959660.734838225,  1955543.76027363), class = c("XY", "POINT", "sfg")), structure(c(951483.685462513,  1939982.16629396), class = c("XY", "POINT", "sfg")), structure(c(951527.767554883,  1940216.55143212), class = c("XY", "POINT", "sfg")), structure(c(958310.673042469,  1954704.68186897), class = c("XY", "POINT", "sfg")), structure(c(950492.05212104,  1951434.68524296), class = c("XY", "POINT", "sfg")), structure(c(959660.734838225,  1955543.76027363), class = c("XY", "POINT", "sfg")), structure(c(959207.877070254,  1944457.42827898), class = c("XY", "POINT", "sfg")), structure(c(960500.020456073,  1955292.64874361), class = c("XY", "POINT", "sfg")), structure(c(959660.734838225,  1955543.76027363), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT",  "sfc"), precision = 0, bbox = structure(c(xmin = 950492.05212104,  ymin = 1939982.16629396, xmax = 960500.020456073, ymax = 1955543.76027363 ), class = "bbox"), crs = structure(list(input = "EPSG:5179",      wkt = "PROJCRS[\"Korea 2000 / Unified CS\",\n    BASEGEOGCRS[\"Korea 2000\",\n        DATUM[\"Geocentric datum of Korea\",\n            ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4737]],\n    CONVERSION[\"Korea Unified Belt\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",38,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",127.5,\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\",1000000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",2000000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"northing (X)\",north,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"easting (Y)\",east,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"unknown\"],\n        AREA[\"Korea, Republic of (South Korea)\"],\n        BBOX[28.6,122.71,40.27,134.28]],\n    ID[\"EPSG\",5179]]"), class = "crs"), n_empty = 0L)), row.names = c(4177L,  15721L, 21365L, 21973L, 24836L, 59359L, 66313L, 70379L, 83277L,  90828L), class = c("sf", "data.frame"), sf_column = "geometry", agr = structure(c(X = NA_integer_,  Y = NA_integer_), .Label = c("constant", "aggregate", "identity" ), class = "factor"))

library(spatstat)

#Convert to ppp format: method 1
pp1 <- as.ppp(pp)
#Warning message:
#In as.ppp.sf(pp) : only first attribute column is used for marks
pp1 <- unmark(pp1)

#Convert to ppp format: method 2
pp2 <- as.ppp(st_coordinates(pp), st_convex_hull(st_union(pp)))
#Warning message:
#data contain duplicated points

identical(pp1, pp2)
[1] FALSE

The first method requires marks and I haven't been able to find out how to turn that off (specifying marks=NULL does not work). I remove marks afterwards but that is a workaround.

The second method is probably wrong, considering the warning, although I based it on this solution. The 2 outputs are not identical.

Is any of these methods correct? Why does method 2 produce duplicated points? Is there a more straightforward method of conversion, with no workarounds?


Solution

  • I'm not sure that the following answer works for every version of spatstat, but I think that if you want to generate an unmarked ppp object, you should run

    # packages
    suppressPackageStartupMessages({
      library(sf)
      library(spatstat)
    })
    
    # Initialise sf object
    pp <- structure(list(X = c(959207.877070254, 959660.734838225, 951483.685462513,  951527.767554883, 958310.673042469, 950492.05212104, 959660.734838225,  959207.877070254, 960500.020456073, 959660.734838225), Y = c(1944457.42827898,  1955543.76027363, 1939982.16629396, 1940216.55143212, 1954704.68186897,  1951434.68524296, 1955543.76027363, 1944457.42827898, 1955292.64874361,  1955543.76027363), geometry = structure(list(structure(c(959207.877070254,  1944457.42827898), class = c("XY", "POINT", "sfg")), structure(c(959660.734838225,  1955543.76027363), class = c("XY", "POINT", "sfg")), structure(c(951483.685462513,  1939982.16629396), class = c("XY", "POINT", "sfg")), structure(c(951527.767554883,  1940216.55143212), class = c("XY", "POINT", "sfg")), structure(c(958310.673042469,  1954704.68186897), class = c("XY", "POINT", "sfg")), structure(c(950492.05212104,  1951434.68524296), class = c("XY", "POINT", "sfg")), structure(c(959660.734838225,  1955543.76027363), class = c("XY", "POINT", "sfg")), structure(c(959207.877070254,  1944457.42827898), class = c("XY", "POINT", "sfg")), structure(c(960500.020456073,  1955292.64874361), class = c("XY", "POINT", "sfg")), structure(c(959660.734838225,  1955543.76027363), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT",  "sfc"), precision = 0, bbox = structure(c(xmin = 950492.05212104,  ymin = 1939982.16629396, xmax = 960500.020456073, ymax = 1955543.76027363 ), class = "bbox"), crs = structure(list(input = "EPSG:5179",      wkt = "PROJCRS[\"Korea 2000 / Unified CS\",\n    BASEGEOGCRS[\"Korea 2000\",\n        DATUM[\"Geocentric datum of Korea\",\n            ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4737]],\n    CONVERSION[\"Korea Unified Belt\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",38,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",127.5,\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\",1000000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",2000000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"northing (X)\",north,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"easting (Y)\",east,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"unknown\"],\n        AREA[\"Korea, Republic of (South Korea)\"],\n        BBOX[28.6,122.71,40.27,134.28]],\n    ID[\"EPSG\",5179]]"), class = "crs"), n_empty = 0L)), row.names = c(4177L,  15721L, 21365L, 21973L, 24836L, 59359L, 66313L, 70379L, 83277L,  90828L), class = c("sf", "data.frame"), sf_column = "geometry", agr = structure(c(X = NA_integer_,  Y = NA_integer_), .Label = c("constant", "aggregate", "identity" ), class = "factor"))
    
    # Convert to ppp format: method 1
    (pp1 <- as.ppp(st_geometry(pp)))
    #> Planar point pattern: 10 points
    #> window: rectangle = [950492.1, 960500] x [1939982.2, 1955543.8] units
    

    Created on 2021-03-01 by the reprex package (v0.3.0)

    The two methods that you present are not equivalent because the as.ppp.sfc function defined in the sf package uses a rectangular owin object:

    # packages
    suppressPackageStartupMessages({
      library(sf)
      library(spatstat)
    })
    
    #Initialise sf object
    pp <- structure(list(X = c(959207.877070254, 959660.734838225, 951483.685462513,  951527.767554883, 958310.673042469, 950492.05212104, 959660.734838225,  959207.877070254, 960500.020456073, 959660.734838225), Y = c(1944457.42827898,  1955543.76027363, 1939982.16629396, 1940216.55143212, 1954704.68186897,  1951434.68524296, 1955543.76027363, 1944457.42827898, 1955292.64874361,  1955543.76027363), geometry = structure(list(structure(c(959207.877070254,  1944457.42827898), class = c("XY", "POINT", "sfg")), structure(c(959660.734838225,  1955543.76027363), class = c("XY", "POINT", "sfg")), structure(c(951483.685462513,  1939982.16629396), class = c("XY", "POINT", "sfg")), structure(c(951527.767554883,  1940216.55143212), class = c("XY", "POINT", "sfg")), structure(c(958310.673042469,  1954704.68186897), class = c("XY", "POINT", "sfg")), structure(c(950492.05212104,  1951434.68524296), class = c("XY", "POINT", "sfg")), structure(c(959660.734838225,  1955543.76027363), class = c("XY", "POINT", "sfg")), structure(c(959207.877070254,  1944457.42827898), class = c("XY", "POINT", "sfg")), structure(c(960500.020456073,  1955292.64874361), class = c("XY", "POINT", "sfg")), structure(c(959660.734838225,  1955543.76027363), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT",  "sfc"), precision = 0, bbox = structure(c(xmin = 950492.05212104,  ymin = 1939982.16629396, xmax = 960500.020456073, ymax = 1955543.76027363 ), class = "bbox"), crs = structure(list(input = "EPSG:5179",      wkt = "PROJCRS[\"Korea 2000 / Unified CS\",\n    BASEGEOGCRS[\"Korea 2000\",\n        DATUM[\"Geocentric datum of Korea\",\n            ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4737]],\n    CONVERSION[\"Korea Unified Belt\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",38,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",127.5,\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\",1000000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",2000000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"northing (X)\",north,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"easting (Y)\",east,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"unknown\"],\n        AREA[\"Korea, Republic of (South Korea)\"],\n        BBOX[28.6,122.71,40.27,134.28]],\n    ID[\"EPSG\",5179]]"), class = "crs"), n_empty = 0L)), row.names = c(4177L,  15721L, 21365L, 21973L, 24836L, 59359L, 66313L, 70379L, 83277L,  90828L), class = c("sf", "data.frame"), sf_column = "geometry", agr = structure(c(X = NA_integer_,  Y = NA_integer_), .Label = c("constant", "aggregate", "identity" ), class = "factor"))
    
    # Convert to ppp format: method 1
    (pp1 <- as.ppp(st_geometry(pp)))
    #> Planar point pattern: 10 points
    #> window: rectangle = [950492.1, 960500] x [1939982.2, 1955543.8] units
    
    # Method 2
    (pp2 <- as.ppp(
      X = st_coordinates(pp), 
      W = as.owin(st_bbox(pp))
    ))
    #> Warning: data contain duplicated points
    #> Planar point pattern: 10 points
    #> window: rectangle = [950492.1, 960500] x [1939982.2, 1955543.8] units
    
    identical(pp1, pp2)
    #> [1] TRUE
    

    Moreover, the first method returns no warning message since the authors set check = FALSE.

    anyDuplicated(pp1)
    #> [1] TRUE
    

    Created on 2021-03-01 by the reprex package (v0.3.0)

    See here for more details.