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.
Best Answer
Calling
st_as_sf
on something that is already ansf
object appears to be a null operation:Sample data frame:
Make into an sf object in 4326 using x and y:
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:oh that looks identical to
s
. If I really want to create ansf
from an existingsf
using different columns and a different CRS then I have to wrap it inas.data.frame
to first remove itssf
-ishness:So note that's now the
a
andb
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:is a no-op -
gw_final
starts as ansf
object (fromst_intersection
) so it doesn't change. At that point anyway it doesn't have columns calledLATITUDE
orLONGITUDE
so I'm not sure what you are expecting it to do anyway. Did you really want tost_transform
it to that CRS?Anyway, with that unchanged, its still in lat-long when you pass it to the
as
function, andspatstat
is only happy with cartesian projected coordinates, which lat-long isn't, so it errors.