I am new to spatstat package and I am trying to make a point patten map with a shape file and some case points. I was able to plot the shape file plot but I just don't know how to make the points on the map plot through the spastat object even I tried points() function or plot() function. What am I missing?

pacman::p_load(sf, sp, spatstat)

# shape file
manchester_ward <- st_read("https://raw.githubusercontent.com/RUMgroup/Spatial-data-in-R/master/rumgroup/data/wards.geojson")
Fallowfield <- manchester_ward %>% dplyr::filter(wd16nm == "Fallowfield")
Fallowfield_pp <- as.owin(Fallowfield)

# points
crime <- structure(list(crime_type = c("Burglary", "Burglary", "Burglary", 
"Burglary", "Burglary", "Burglary"), geometry = structure(list(
    structure(c(-2.242511, 53.44558), class = c("XY", "POINT", 
    "sfg")), structure(c(-2.237017, 53.445735), class = c("XY", 
    "POINT", "sfg")), structure(c(-2.2334, 53.442408), class = c("XY", 
    "POINT", "sfg")), structure(c(-2.242048, 53.448871), class = c("XY", 
    "POINT", "sfg")), structure(c(-2.242048, 53.448871), class = c("XY", 
    "POINT", "sfg")), structure(c(-2.243242, 53.444339), class = c("XY", 
    "POINT", "sfg"))), class = c("sfc_POINT", "sfc"), precision = 0, bbox = structure(c(xmin = -2.243242, 
ymin = 53.442408, xmax = -2.2334, ymax = 53.448871), class = "bbox"), crs = structure(list(
    input = "EPSG:4326", wkt = "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"Horizontal component of 3D system.\"],\n        AREA[\"World.\"],\n        BBOX[-90,-180,90,180]],\n    ID[\"EPSG\",4326]]"), class = "crs"), n_empty = 0L)), sf_column = "geometry", agr = structure(c(crime_type = 1L), .Label = c("constant", 
"aggregate", "identity"), class = "factor"), row.names = c(NA, 
6L), class = c("sf", "data.table", "data.frame"))

crime_pt <- matrix(unlist(crime$geometry), ncol = 2, byrow = T)
k1 <- ppp(x = crime_pt[,1], y = crime_pt[,2],
          window = Fallowfield_pp, check = T)

  • What's windows in the ppp call at the end? Also you should never load tidyverse since it pulls in so much - all you are using from that is filter which you could get better with dplyr::filter (since then it couldn't clash with base::filter). – Spacedman Mar 19 at 12:29
  • Thanks @Spacedman. The windows was supposed to be "Fallowfield_pp". I edited out the tidyverse and others. I got a warning message: Warning message: In plot.ppp(k1) : 6 illegal points also plotted and an almost blank page in the plot panel. – Gnin Mar 19 at 12:35

Different coordinate systems. Fallowfield is this:

> Fallowfield_pp 
window: polygonal boundary
enclosing rectangle: [382951.5, 385869.8] x [393616.3, 394988.8] units

and your crime is:

> crime
Simple feature collection with 6 features and 1 field
Attribute-geometry relationship: 1 constant, 0 aggregate, 0 identity
geometry type:  POINT
dimension:      XY
bbox:           xmin: -2.243242 ymin: 53.44241 xmax: -2.2334 ymax: 53.44887
geographic CRS: WGS 84
  crime_type                   geometry
1   Burglary POINT (-2.242511 53.44558)
2   Burglary POINT (-2.237017 53.44573)
3   Burglary   POINT (-2.2334 53.44241)
4   Burglary POINT (-2.242048 53.44887)

notice the very different coordinate numbers. Fallowfield is British National Grid and crime is WGS84 lat-long.

Let's convert the crime data to the same coordinate system as the Fallowfield boundary:

 crime = st_transform(crime, st_crs(Fallowfield))

and then go through your ppp construction again and you get:

enter image description here

with a duplicate point warning, because point pattern statistics is based on points being generated on the continuum, but that's not important right now I guess...

  • Thanks and I didn't know they're in different systems! And I didn't recognise you earlier - big fan of your spatial course! I guess I should have given you more points! – Gnin Mar 19 at 13:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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