Quadrats in Spatstat

Question

I am now doing exploratory analysis and the objective is to plot a quadrat map, determine if there exists complete spatial randomness (visually and using chi-square test), etc. However, I am having trouble plotting a quadrat map. I previously asked this question that led me to reproject my data.

Here is my code:

library(rgdal) #Brings Spatial Data in R
library(spatstat) # Spatial Statistics
library(maptools)
library(raster)

# Load nyc zip code boundary polygon shapefile 
s <- readOGR("/Users/my_name/Documents/fproject/zip","zip")
nyc <- as(s,"owin")

### OGR data source with driver: ESRI Shapefile 
#Source: "/Users/my_name/Documents/project/zip", layer: "zip"
#with 263 features

# Load nyc arrests point feature shapefile
s <- readOGR("/Users/my_name/Documents/project/nycarrests/","geo1")

### OGR data source with driver: ESRI Shapefile 
#Source: "/Users/my_name/Documents/project/nycarrests", layer: "geo1"
#with 103376 features
#It has 19 fields

#Converting the dataset into a point pattern
utm <- "+proj=utm +zone=32 +datum=WGS84"
x <- spTransform(s,utm)
arrests <- as(x,"ppp")

marks(arrests) <- NULL
Window(arrests) <- nyc 
plot(arrests, main="2020", cols="dark green",pch=20)

### Quadrat density 2020 
Q <- quadratcount(arrests, nx= 3, ny=3)
plot(arrests, pch=20, cols="grey70", main="2020")  # Plot points
plot(Q, main="2020", add=TRUE)  # Add quadrat grid

Here, my data points are constrained by NYC zip code boundaries, but they do not show up on my map. What can I do to fix this?

After you create each object, I suggest that you print the object (by just typing its name), and look at the output to see whether it is what you expected. For example, after creating `arrests`, print it to see how many points have been included, and read the description of the spatial domain to see if it is what you expected. You will probably find that you've forgotten to spatially transform some of the data. — Adrian Baddeley, Dec 19 '20 at 14:34

score 0 · Answer 1 · answered Dec 15 '20 at 22:14

You are not providing much information. For example you could show (print) s (each time), or crs(s) (both with the raster package loaded)

Either way, it is almost certainly because the coordinate reference systems you are using are not the same. You change one to "UTM zone 21", but you do not change the others. They are probably different. Also, zone 21 is not appropriate for New York.

Quadrats in Spatstat

1 Answers1