I have a question similar to this one from 2014, which was answered but the datasets are no longer available and our original data structures differ. (I'm in crunch time and stumped, so if you're able to respond quickly I would greatly appreciate it!!)
Goal: use the type of bedrock as a covariate in a Point Process Model (ppm) in spatstat with mine locations in Connecticut.
Data: the files are available from this Dropbox folder. The rock data and CT poly outline comes from UConn Magic Library, and the mine data comes from the USGS Mineral Resources Data System.
Approach: I loaded some relevant packages and read in the shapefiles (and converted coords to match CT's system), and used the CT polygon as an owin object.
library(rgdal)
library(splancs)
library(spatstat)
library(sp)
library(raster)
library(geostatsp)
#read in shapefiles
ct <-readOGR(".","CONNECTICUT_STATE_POLY")
mrds <-readOGR(".","mrds-2017-02-20-23-30-58")
rock<-readOGR(".","bedrockpolyct_37800_0000_2000_s50_ctgnhs_1_shp_wgs84")
#convert mrds and rock to ct's coord system
tempcrs<-ct@proj4string
mrds<-spTransform(mrds,tempcrs)
rock<-spTransform(rock,tempcrs)
#turn ct shapefile into owin, call it w for window
w <-as.owin(ct)
#subset mrds data to just CT mines
mrdsCT <-subset(mrds,mrds@data$state=="Connecticut")
#ppm can't handle marked data yet, so need to unmark()
#create ppp object for mrds data, set window to w
mrdsCT.ppp <-as.ppp(mrdsCT)
Window(mrdsCT.ppp)<-w
From "Modelling Spatial Point Patterns in R" by Baddeley & Turner (page 39): Unfortunately a pixel image in spatstat cannot have categorical (factor) values, because R refuses to create a factor-valued matrix. In order to represent a categorical variate as a pixel image, the categorical values should be encoded as integers (for efficiency’s sake) and assigned to an integer-valued pixel image. Then the model formula should invoke the factor command on this image. For example if fim is an image with integer values which represent levels of a factor, then:
ppm(X, ˜factor(f), Poisson(), covariates=list(f=fim))
There are several different types of rock classification included in the shapefile. I'm interested in LITHO1, which is a factor with 27 levels. It's the sixth attribute.
litho1<-rock[,6]
My (limited but researched) understanding is that I need to convert the shapefile to a raster, and later convert it to an image in order to be used in ppm. I created a mask from ct, and used that.
ctmask<-raster(ct, resolution=2000)
ctmask[!is.na(ctmask)] <- 0
litho1rast<-rasterize(litho1,ctmask)
After this point, I've tried several approaches and haven't had success just yet. I've attempted to follow the approaches laid out in the question linked, as well as search in documentation for relevant examples to adopt (factor, ratify, levels). Unlike the prior question, my data was already a factor, so it wasn't clear why I should apply the factor function to it.
Looking at litho1rast, the @data@attributes dataframe contains the following. If I plot it, it just plots the ID; levelplot function does plot LITHO1. When I would apply the factor functions, the ID would be retained but not LITHO1.
$ ID : int [1:1891] 1 2 3 4 5 6 7 8 9 10 ...
$ LITHO1: Factor w/ 27 levels "amphibolite",..: 23 16 23 16 23 16 24 23 16 24 ...
The ppm model would need an object class im, so I converted the raster to the im. I tried two ways. I can make ppm execute...but it treats every point as a factor rather than the 27 levels (with either litho1.im or litho1.im2) ...
litho1.im<-as.im(litho1rast)
litho1.im2<-as.im.RasterLayer(litho1rast)
model1=ppm(unmark(mrdsCT.ppp) ~ factor(COV1), covariates=list(COV1=litho1.im))
model1
So, I'm not quite sure where to go from here. It seems like I need to pass an argument to the as.im so that it knows to retain the LITHO1 not the ID. Clever ideas or leads to pertinent functions or approaches much appreciated!
