Creating a point pattern (ppp) error using SpatStat in R: "Only projected coordinates may be converted to spatstat class objects"

Question

I have a point dataset, and I am taked with filtering out duplicate coordinates, removing all points that are outside of a specified county shapefile, and then creating a ppp points object. However, at the end of this operation, I get the following error:

Error in as.ppp.SpatialPointsDataFrame(from) : 
  Only projected coordinates may be converted to spatstat class objects

How can I fix this error?

The two datasets you will need to reproduce this error can be found here (a .csv file and a shape file).

Here is what I have tried:

#Load packages and datasets
library(pacman)
p_load(spatstat,
       dplyr,
       maptools,
       raster,
       sf,
       sp,
       ggplot2)

gw <- read.csv("RileyCNTYGWwells.csv", stringsAsFactors = FALSE)
KS_counties <- st_read("KS_counties.shp")

#Select relevant columns
gw_sp <- gw %>%
  dplyr::select(LONGITUDE, LATITUDE, WELL_USE, WELL_DEPTH, EST_YIELD) %>% na.omit(gw_sp)

#Project to WGS84
gw_cor <- st_as_sf(gw_sp, coords=c("LONGITUDE","LATITUDE"),
                crs = st_crs(4326))

#Remove duplicate coordinates with dplyr's "distinct" function
gw_sp <- gw_cor %>% 
  distinct()

crs(gw_sp)

#Filter KS shapefile to desired county
riley <- KS_counties %>% 
  filter(name == "Riley")

#Transform to WGS84, same as the gw_sp object
riley <- st_transform(riley, 4326)

#Filter data that is outside of the county line. It may be necessary to highlight both of these lines and run them simultaneously to get this code to work.
plot(riley$geometry)
plot(gw_sp$geometry, add = T)

#Get points within the county
gw_final <- st_intersection(riley, gw_sp)

plot(riley$geometry, main = "Final filtered points")
plot(gw_final$geometry, add = T)

#Convert back to ppp object
gw_final <- st_as_sf(gw_final, coords=c("LONGITUDE","LATITUDE"),
                crs = crs('+proj=lcc +lat_1=38.56666666666667 +lat_2=37.26666666666667 +lat_0=36.66666666666666 +lon_0=-98.5 +x_0=400000 +y_0=400000 +ellps=GRS80 +datum=NAD83 +units=m +no_defs'))

gw_ppp <- as(gw_final, "Spatial")
gw_ppp <- as(gw_ppp, "ppp")

The error mentioned above occurs after the last line is run.

How can I fix this error?

I thought I was projecting it when I ran the st_as_sf() command above.

Answer 1

Calling st_as_sf on something that is already an sf object appears to be a null operation:

Sample data frame:

> d = data.frame(x=1:5, y=1:5, a=5:1, b=5:1)
> d
  x y a b
1 1 1 5 5
2 2 2 4 4
3 3 3 3 3
4 4 4 2 2
5 5 5 1 1

Make into an sf object in 4326 using x and y:

> s = st_as_sf(d, coords=1:2, crs=4326)
> s
Simple feature collection with 5 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 1 ymin: 1 xmax: 5 ymax: 5
Geodetic CRS:  WGS 84
  a b    geometry
1 5 5 POINT (1 1)
2 4 4 POINT (2 2)
3 3 3 POINT (3 3)
4 2 2 POINT (4 4)
5 1 1 POINT (5 5)

okay good looks like 5 points in WGS84 epsg:4326. Now let's try making another sf from that using the other two columns and a different coord system:

> s2 = st_as_sf(s, coords=c("a","b"), crs=27700)
> s2
Simple feature collection with 5 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 1 ymin: 1 xmax: 5 ymax: 5
Geodetic CRS:  WGS 84
  a b    geometry
1 5 5 POINT (1 1)
2 4 4 POINT (2 2)
3 3 3 POINT (3 3)
4 2 2 POINT (4 4)
5 1 1 POINT (5 5)
> 

oh that looks identical to s. If I really want to create an sf from an existing sf using different columns and a different CRS then I have to wrap it in as.data.frame to first remove its sf-ishness:

> s2 = st_as_sf(as.data.frame(s), coords=c("a","b"), crs=27700)
> s2
Simple feature collection with 5 features and 0 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 1 ymin: 1 xmax: 5 ymax: 5
Projected CRS: OSGB 1936 / British National Grid
     geometry
1 POINT (5 5)
2 POINT (4 4)
3 POINT (3 3)
4 POINT (2 2)
5 POINT (1 1)

So note that's now the a and b columns in OSGB BNG coords. Note also the coords are the same numbers, this is not doing a transform, its just assigning the projection code to those numbers. So your line:

gw_final <- st_as_sf(gw_final, coords=c("LONGITUDE","LATITUDE"),
                crs = crs('...')

is a no-op - gw_final starts as an sf object (from st_intersection) so it doesn't change. At that point anyway it doesn't have columns called LATITUDE or LONGITUDE so I'm not sure what you are expecting it to do anyway. Did you really want to st_transform it to that CRS?

Anyway, with that unchanged, its still in lat-long when you pass it to the as function, and spatstat is only happy with cartesian projected coordinates, which lat-long isn't, so it errors.