1

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.

1 Answer 1

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.

1
  • Thanks for your help! Really appreciate it. Apr 19, 2022 at 22:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.