R包spatstat：如何使用点过程模型协变量作为从形状文件开始的因素
R package spatstat: How to use point process model covariate as factor starting with shape file
我有一个类似于 this one from 2014, 的问题已得到解答，但数据集不再可用且我们的原始数据结构不同。 （我在紧要关头被难住了，所以如果你能迅速回复我将不胜感激！！）
目标：在 spatstat 中的点过程模型 (ppm) 中使用基岩类型作为协变量，矿区位于康涅狄格州。
数据：文件可从Dropbox folder获得。岩石数据和CT多边形轮廓来自UConn Magic Library，矿山数据来自USGS Mineral Resources Data System
方法：
我加载了一些相关的包并读入了 shapefile（并转换坐标以匹配 CT 的系统），并将 CT 多边形用作 owin 对象。
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
来自 Baddeley & Turner 的 "Modelling Spatial Point Patterns in R"（第 39 页）：
不幸的是，spatstat 中的像素图像不能具有分类（因子）值，因为 R 拒绝创建因子值矩阵。为了将分类变量表示为像素图像，分类值应编码为整数（出于效率考虑）并分配给整数值像素图像。然后模型公式应该在这个图像上调用 factor 命令。例如，如果 fim 是一个具有表示因子水平的整数值的图像，则：
ppm(X, ˜factor(f), Poisson(), covariates=list(f=fim))
shapefile 中包含几种不同类型的岩石 classification。我对 LITHO1 很感兴趣，它是一个具有 27 个级别的因子。这是第六个属性。
litho1<-rock[,6]
我（有限但经过研究）的理解是我需要将 shapefile 转换为光栅，然后再将其转换为图像以便在 ppm 中使用。我从 ct 创建了一个面具，并使用了它。
ctmask<-raster(ct, resolution=2000)
ctmask[!is.na(ctmask)] <- 0
litho1rast<-rasterize(litho1,ctmask)
在这一点之后，我尝试了几种方法，但还没有成功。我试图遵循链接问题中列出的方法，并在文档中搜索要采用的相关示例（因子、批准、级别）。与前面的问题不同，我的数据已经是一个因素，所以不清楚为什么我应该对它应用因素函数。
查看 litho1rast，@data@attributes 数据框包含以下内容。如果我绘制它，它只会绘制 ID； levelplot 函数绘制 LITHO1。当我应用因子函数时，ID 会被保留，但 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 ...
ppm 模型需要一个对象 class im，因此我将栅格转换为 im。我尝试了两种方法。我可以让 ppm 执行...但它将每个点都视为一个因素，而不是 27 个级别（使用 litho1.im 或 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
所以，我不太确定从这里到哪里去。似乎我需要向 as.im 传递一个参数，以便它知道保留 LITHO1 而不是 ID。 非常感谢聪明的想法或导致相关功能或方法！
查看您的代码，您似乎没有向 rasterize 提供 field 参数。
来自rasterize帮助：
fieldnumeric or character. The value(s) to be transferred. This can
be a single number, or a vector of numbers that has the same length as
the number of spatial features (points, lines, polygons). If x is a
Spatial*DataFrame, this can be the column name of the variable to be
transferred. If missing, the attribute index is used (i.e. numbers
from 1 to the number of features). You can also provide a vector with
the same length as the number of spatial features, or a matrix where
the number of rows matches the number of spatial features
在这一行：
litho1rast<-rasterize(litho1,ctmask)
您可能必须指定 litho 对象的哪一列用于光栅化。类似于：
litho1rast<-rasterize(litho1,ctmask, field = "LITHO1")
引用的 Baddeley & Turner 声明不再正确 --- 该引用来自一套非常古老的车间笔记。
classim 的像素图像可以有因子值（自 2011 年起）。如果 Z 是整数值像素图像（class im），您可以通过设置 levels(Z) <- lev where lev 使其成为因子值图像是可能值的标签的特征向量。
你应该不需要使用rasterize：应该可以使用as.im将rock[,6]直接转换成像素图像（加载maptools包后） .
请参阅 Baddeley、Rubak 和 Turner 的书（空间点模式：R 的方法和应用，CRC 出版社，2016 年）以获得完整解释。
我有一个类似于 this one from 2014, 的问题已得到解答，但数据集不再可用且我们的原始数据结构不同。 （我在紧要关头被难住了，所以如果你能迅速回复我将不胜感激！！）
目标：在 spatstat 中的点过程模型 (ppm) 中使用基岩类型作为协变量，矿区位于康涅狄格州。
数据：文件可从Dropbox folder获得。岩石数据和CT多边形轮廓来自UConn Magic Library，矿山数据来自USGS Mineral Resources Data System
方法： 我加载了一些相关的包并读入了 shapefile（并转换坐标以匹配 CT 的系统），并将 CT 多边形用作 owin 对象。
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
来自 Baddeley & Turner 的 "Modelling Spatial Point Patterns in R"（第 39 页）： 不幸的是，spatstat 中的像素图像不能具有分类（因子）值，因为 R 拒绝创建因子值矩阵。为了将分类变量表示为像素图像，分类值应编码为整数（出于效率考虑）并分配给整数值像素图像。然后模型公式应该在这个图像上调用 factor 命令。例如，如果 fim 是一个具有表示因子水平的整数值的图像，则：
ppm(X, ˜factor(f), Poisson(), covariates=list(f=fim))
shapefile 中包含几种不同类型的岩石 classification。我对 LITHO1 很感兴趣，它是一个具有 27 个级别的因子。这是第六个属性。
litho1<-rock[,6]
我（有限但经过研究）的理解是我需要将 shapefile 转换为光栅，然后再将其转换为图像以便在 ppm 中使用。我从 ct 创建了一个面具，并使用了它。
ctmask<-raster(ct, resolution=2000)
ctmask[!is.na(ctmask)] <- 0
litho1rast<-rasterize(litho1,ctmask)
在这一点之后，我尝试了几种方法，但还没有成功。我试图遵循链接问题中列出的方法，并在文档中搜索要采用的相关示例（因子、批准、级别）。与前面的问题不同，我的数据已经是一个因素，所以不清楚为什么我应该对它应用因素函数。
查看 litho1rast，@data@attributes 数据框包含以下内容。如果我绘制它，它只会绘制 ID； levelplot 函数绘制 LITHO1。当我应用因子函数时，ID 会被保留，但 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 ...
ppm 模型需要一个对象 class im，因此我将栅格转换为 im。我尝试了两种方法。我可以让 ppm 执行...但它将每个点都视为一个因素，而不是 27 个级别（使用 litho1.im 或 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
所以，我不太确定从这里到哪里去。似乎我需要向 as.im 传递一个参数，以便它知道保留 LITHO1 而不是 ID。 非常感谢聪明的想法或导致相关功能或方法！
查看您的代码，您似乎没有向 rasterize 提供 field 参数。
来自rasterize帮助：
fieldnumeric or character. The value(s) to be transferred. This can be a single number, or a vector of numbers that has the same length as the number of spatial features (points, lines, polygons). If x is a Spatial*DataFrame, this can be the column name of the variable to be transferred. If missing, the attribute index is used (i.e. numbers from 1 to the number of features). You can also provide a vector with the same length as the number of spatial features, or a matrix where the number of rows matches the number of spatial features
在这一行：
litho1rast<-rasterize(litho1,ctmask)
您可能必须指定 litho 对象的哪一列用于光栅化。类似于：
litho1rast<-rasterize(litho1,ctmask, field = "LITHO1")
引用的 Baddeley & Turner 声明不再正确 --- 该引用来自一套非常古老的车间笔记。
classim 的像素图像可以有因子值（自 2011 年起）。如果 Z 是整数值像素图像（class im），您可以通过设置 levels(Z) <- lev where lev 使其成为因子值图像是可能值的标签的特征向量。
你应该不需要使用rasterize：应该可以使用as.im将rock[,6]直接转换成像素图像（加载maptools包后） .
请参阅 Baddeley、Rubak 和 Turner 的书（空间点模式：R 的方法和应用，CRC 出版社，2016 年）以获得完整解释。