CODEWITHSUNDEEP
rcoordinate-systemsr-sfspatstat
SavedByJESUS
SavedByJESUS

Reputation: 3314

Project longitude and latitude to planar coordinates system

Disclaimer: I am currently using the ppp branch version of the {sf} package, because new features for converting objects between {sf} and {spatstat} are available in it (see https://github.com/r-spatial/sf/issues/1233). For it to work properly, I had to manually delete the {sf} package from my hard drive and then reinstall it from Github. I am also using the development version of {spatstat} for no particular reason.

# install.packages("remotes")
install_github("r-spatial/sf@ppp")
install_github("spatstat/spatstat")

I have two geospatial objects: area_one, which is the union of the polygons of several counties in Texas and vz, which is the point locations of several stores in Texas and they are both objects of the sf family. vz was created using longitude and latitude coordinates scraped from the internet. My goal is to create a ppp object with the locations in vz as the points and the polygon in area_one as the window. The issue is that I cannot find the correct coordinate reference system (CRS), for my points to lie inside the polygon. I get an error telling me that the points lie outside the window. Here are the two files to make the code below reproducible:

area_one: download here

vz: download here

# Load packages

library(sf) # Development version in the ppp branch
library(spatstat) # Development version in the master branch
library(tmap)
library(here)

# Read the geospatial data (CRS = 4326)

area_one <- st_read(dsn = here("area_one/area_one.shp"), layer = "area_one")
vz <- st_read(dsn = here("vz/vz.shp"), layer = "vz")

# Plot a quick map

tm_shape(area_one) +
  tm_borders() +
  tm_shape(vz) +
  tm_dots(col = "red", size = 0.1)

enter image description here

# Create a planar point pattern

vz_lonlat <- st_coordinates(vz)
area_one_flat <- st_transform(area_one, crs = 6345)

p <- ppp(x = vz_lonlat[, "X"], y = vz_lonlat[, "Y"], window = as.owin(area_one_flat)) # Note that the reason why this line of code does not throw an error (especially despite the window argument being an sf object) is because of the version of the {sf} package I am using

Warning message:
49 points were rejected as lying outside the specified window 

plot(p)

enter image description here

Upvotes: 1

Views: 775

Answers (1)

Ege Rubak
Ege Rubak

Reputation: 4507

As @spacedman points out you should first transform vz to the same coordinate system as the observation region. I guess you could do something like (untested):

vz_flat <- st_coordinates(st_transform(vz, crs = 6345))
area_one_flat <- st_transform(area_one, crs = 6345)
p <- ppp(x = vz_flat[, "X"], y = vz_flat[, "Y"], window = as.owin(area_one_flat))

Upvotes: 1

Related Questions