创建具有相同范围的多个密度图

Question

我正在尝试使用tmaptools 包中的smooth_map 函数从点创建密度图（栅格）。我想在我的数据框 (df) 中为变量 (var) 的每个因子水平创建一个栅格。所有栅格都应具有相同的范围，由 shapefile (study_area) 给出。我想堆叠它们，然后使用levelplot 绘制它们。

但是，该函数返回的所有栅格的范围略有不同。我不知道如何创建具有相同范围的密度图。这是代码：

library(tmaptools) 
library(rgdal)
library(sf)
library(raster)

##reproducible example of dataframe
df<-setNames(data.frame(matrix(ncol = 3, nrow = 7)), c("longitude", "latitude", "var"))

df$longitude<-c(-53.30002, -50.47749, -55.85561, -57.88447, -55.83864, -58.84610, -57.49215)
df$latitude<-c(-13.037530,  -13.023480, -13.416200, -13.659120, -12.758670, -13.114460, -14.622520)
df$var<-c("A", "A", "A", "A", "B", "B", "B")

## convert var to factor
df$var<-as.factor(df$var)

##read shapefile of Mato Grosso, Brazil 
## shapefile can be downloaded from e.g. https://geo.nyu.edu/catalog/stanford-vw409fv3488
study_area<-readOGR("shape_MT.shp", layer="shape_MT") #will load the shapefile 

#create empty stack
s <- stack()

#loop over factor levels
for (i in levels(df$var)) { 
  a<-subset(df, df$var==i)
  my.sf.point <- st_as_sf(x = a, coords = c("longitude", "latitude"),crs = "+proj=longlat  +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
  density_map <- smooth_map(my.sf.point, cover = study_area) 
  r<-density_map$raster
  print(extent(r))
  s <- stack(s,r)
}

我收到以下错误消息：

compareRaster(x) 中的错误：不同程度

关于如何解决这个问题的任何想法？

Answer 1

也许这对你有用：

## convert var to factor
df$var<-as.factor(df$var)

## read shapefile of Mato Grosso, Brazil 
## shapefile can be downloaded from e.g. https://geo.nyu.edu/catalog/stanford-vw409fv3488
study_area<-readOGR("....", layer="...") 

它也可以与您自己的 crs 一起使用，但也许 ESRI:102033 South_America_Albers_Equal_Area_Conic 在这里看起来是一个更好的选择（或类似的东西）：

crs.laea <- CRS("+proj=aea +lat_1=-5 +lat_2=-42 +lat_0=-32 +lon_0=-60 +x_0=0 +y_0=0 +ellps=aust_SA +units=m +no_defs")
study_area <- spTransform(study_area, CRS = crs.laea)

#create empty stack
s <- list()

#loop over factor levels
for (i in levels(df$var)) { 
  a<-subset(df, df$var==i)
  my.sf.point <- st_as_sf(x = a, coords = c("longitude", "latitude"), crs = "+proj=longlat +ellps=aust_SA +no_defs")

  #reprojecting to ESRI:102033 South_America_Albers_Equal_Area_Conic
  my.sf.point <- st_transform(my.sf.point, 102033)
  density_map <- smooth_map(my.sf.point, cover = study_area, nrow = 514, ncol = 510) 
  r<-density_map$raster

  extent(r) <- c(-175916.2, 1048374, 1610251, 2845008) # extent of the study area
  plot(r)

  s[[i]] <- r
}

s[["A"]] <- resample(s[["A"]],s[["B"]],method='bilinear')

stack(s)