如何使用ggplot2中的shapefiles从第一个图中保留scale_fill_color

Question

我需要使用 ggplot2 中的 shapefile 来进行scale_fill_manual 的帮助。我尝试了很多事情，终于发布了，希望有人能给我提示。我基本上正在绘制一个 shapefile 并使用 scale_fill_manual 用自定义颜色对其进行可视化，然后在其上覆盖一些点，但是当我尝试在图例中包含我的新点时，我的原始颜色在那里，但值都弄乱了。绘制形状文件的上部工作正常，但覆盖新点时的底部是我需要帮助的地方。我内嵌了一些评论。如下：shapefile的下载路径为： https://login.filesanywhere.com/fs/v.aspx ?v=8c6a66875b6574bbaa68

library(tidyverse)
library(rgdal)
library(maptools)
library(plyr)
library(sp)
library(geosphere)
library(data.table)
library(rgeos)

wolves.map <- readOGR(dsn=".", layer="PNW_wolf_habitat_grid")
 message(proj4string(wolves.map)) # it is in Albers Equal Area projection.

 #Select presence/abscense only (1 and 0)
wolfsub <- wolves.map[!wolves.map$WOLVES_99 %in% 2,]
wolfsub$MAJOR_LC <-as.numeric(as.character(wolfsub$MAJOR_LC))

# Add columns to the wolfsub dataset. 42 = Forest, 51 = Shrub, > 81 = Agriculture
wolfsub$Forest<-ifelse(wolfsub$MAJOR_LC==42,1,0)
wolfsub$Shrub<-ifelse(wolfsub$MAJOR_LC==51,1,0)
wolfsub$Agriculture <- ifelse(wolfsub$MAJOR_LC > 81,1,0)

# create the model
mod1<-glm(WOLVES_99 ~ RD_DENSITY + Forest + Shrub + Agriculture,family = binomial,data = wolfsub)
summary(mod1)

#fitted(mod1)
wolfsub$WOLVES_99pred <- fitted(mod1) # add the predicted values to wolfsub

# Convert the wolves.map shapefile to data.frame
wolves.mapDF <- as.data.frame(wolves.map)

#fortify wolves.map to be used with ggplot2
wolves.ds <- fortify(wolves.map,region="GRID2_ID")

 # Rename the 'GRID2'_ID to 'id' to to be able to merge with the shapefile wolves.map
 wolves.mapDF <- rename(wolves.mapDF,c(GRID2_ID="id"))

# merge the shapefile wolves.ds and wolves.mapDF dataframe to be able to use the wolves.mapDF variables with ggplot2
wolves.ggmap <- merge(wolves.ds, wolves.mapDF, by = "id", all = TRUE)
wolves.ggmap <- wolves.ggmap[order(wolves.ggmap$order), ]

 wolves.ggmap$MAJOR_LC <-as.numeric(as.character(wolves.ggmap$MAJOR_LC))

 ### Now do the whole data set
 # 42 = Forest, 51 = Shrub, > 81 = Agriculture
wolves.ggmap$Forest<-ifelse(wolves.ggmap$MAJOR_LC==42,1,0)
wolves.ggmap$Shrub<-ifelse(wolves.ggmap$MAJOR_LC==51,1,0)
wolves.ggmap$Agriculture<-ifelse(wolves.ggmap$MAJOR_LC>81,1,0)

# Predict probabilities for the whole dataset
wolves.ggmap$PredictedSuit <- predict(mod1,newdata=wolves.ggmap,type='response')

#Make PredictedSuit a factor
wolves.ggmap$DiscretePred <- cut(wolves.ggmap$PredictedSuit,breaks=c(0,0.29,0.40,0.45,0.6,0.69),dig.lab = 2,include.lowest=TRUE)

 #plot and display a legend with the new cuts
Palette1 <- c('grey80','orange','yellow','green','green3','blue')
wolves.pred3 <- ggplot(wolves.ggmap,aes(long,lat,group=group)) + theme_bw() + theme_void() +  
geom_polygon(aes(fill=DiscretePred), colour = alpha("white", 1/2), size = 0.2) + theme(legend.position = c(0.14, 0.16)) +
scale_fill_manual(values=Palette1)  + guides(fill=guide_legend(ncol=2,"Predicted\n Suitability\n > 0.45"))  
wolves.pred3 

我得到以下图表（好）：

上述所有代码都按预期工作。我遇到的问题在下面。下面的代码可以很好地覆盖上面同一形状文件的子集中的点。但是，当我尝试将新点添加到图例时，我丢失了scale_fill_manual颜色。

#Extract wolves from 2001 first and overlay them on map
wolfsub_01 <- wolves.map[wolves.map$WOLVES_01 %in% 1,]
wolfsub_01$MAJOR_LC <-as.numeric(as.character(wolfsub_01$MAJOR_LC))

#Get centroids to overlay on existing plot
test <- gCentroid(wolfsub_01, byid = TRUE)

#Convert to dataframe to be used with ggplot2
wolf <- as.data.frame(wolfsub_01) 
test <- as.data.frame(test)  
 wolves_test <- cbind(wolf,test)

#Overlay on existing plot
wolves.pred3 +
geom_point(data=wolves_test,aes(x,y,group=NULL,fill='2001 wolves'),color='blue') 

如果我尝试在我的图例中包含“2001 年狼”，我的颜色将保持正确的顺序。然而，我的传奇值变得一团糟。我尝试用不同的调色板重新排列它们，但这只会使情况变得更糟，因为颜色和标签与相应的颜色不对齐。我还需要帮助从图例中删除点。如何将我的颜色恢复到上面在原始绘图上使用的原始 Palette1？可能是一件简单的事情，但我花了很多时间尝试但无法弄清楚。预先感谢。

我明白这个情节。请注意，所有值都已结束。我需要这些值的顺序与第一个图相同。

编辑：这就是我的情节在幕后显示的内容。第一个图具有以下颜色顺序：

> g <- ggplot_build(wolves.pred3)
> unique(g$data[[1]]["fill"])
      fill
1   grey80
9   orange
115 yellow
241 green3
271  green

我的第二个图的颜色顺序与第一个图不同。我想知道如何使第二个颜色顺序与第一个颜色顺序匹配。

> g <- ggplot_build(a)
    > unique(g$data[[1]]["fill"])
          fill
    1   green3
    9   grey80
    115 orange
    241  green
    271 yellow
    > 

Answer 1

这是我为你尝试过的。我浏览了您的所有代码，并给我的印象是您使数据处理变得复杂。我曾经使用这种sp方法并编写像您一样的代码。我认为这种方法使您在某个地方“扭曲”数据处理（例如，当您使用 merge() 时）。在这里，我以另一种方式编写了您的代码，以实现预期的结果。我在下面的脚本中留下了解释。在我看来，要点是避免一些棘手的数据操作。一种方法是使用sf和tidyverse包。我希望这能帮到您。

library(sf)
library(dplyr)
library(ggplot2)
library(rgeos)

# You can use the sf package to read a shapefile.
wolves.map <- st_read(dsn = ".", layer = "PNW_wolf_habitat_grid")

# Step 1
# Sub data
# Select presence/abscense only (1 and 0)
# You used base R to write your script. The sp class objects do not accept
# tidyverse ways. But sf objects can take tidyverse ways, which makes your life much easier.

wolfsub <- filter(wolves.map, WOLVES_99 != 2) %>% 
            mutate(Forest = if_else(MAJOR_LC == 42, 1, 0),
                   Shrub = if_else(MAJOR_LC == 51, 1, 0),
                   Agriculture = if_else(MAJOR_LC > 81, 1, 0))

# Create the model
mod1 <- glm(WOLVES_99 ~ RD_DENSITY + Forest + Shrub + Agriculture, family = binomial, data = wolfsub)
summary(mod1)

# Fitted(mod1)
wolfsub$WOLVES_99pred <- fitted(mod1) # add the predicted values to wolfsub


# Step 2: Whole data
# Here I can avoid creating a new data frame for ggplot2. I saw that you worked
# to arrange a new data frame with all numbers. But that is not necessary any more.

wolves.map %>% 
mutate(Forest = if_else(MAJOR_LC == 42, 1, 0),
       Shrub = if_else(MAJOR_LC == 51, 1, 0),
       Agriculture = if_else(MAJOR_LC > 81, 1, 0)) -> wolves.map 

wolves.map$PredictedSuit <- predict(mod1,newdata = wolves.map,type = 'response')


mutate(wolves.map,
       DiscretePred = cut(PredictedSuit,
                          breaks = c(0,0.29,0.40,0.45,0.6,0.69),
                          dig.lab = 2,include.lowest = TRUE)) -> out

# Plot and display a legend with the new cuts
Palette1 <- c('grey80','orange','yellow','green','green3','blue')

ggplot() +
geom_sf(data = out, aes(fill = DiscretePred),
        colour = alpha("white", 1/2), size = 0.2) +
scale_fill_manual(values = Palette1)  +
theme_bw() +
theme_void() + 
theme(legend.position = c(0.14, 0.16)) +
guides(fill = guide_legend(ncol = 2,"Predicted\n Suitability\n > 0.45")) -> g

# Step 3
# Extract wolves from 2001 first and overlay them on map

wolfsub_01 <- filter(wolves.map, WOLVES_01 == 1)

# Get centroids to overlay on existing plot. I used st_centroid() instead of Gcentroid(). 
# Then, I added long and lat to the original data frame, `wolfsub_01`.
# I also added a new column for color.

test <- bind_cols(wolfsub_01,
                  as.data.frame(st_coordinates(st_centroid(wolfsub_01)))) %>%
        mutate(color = "blue")

# Finally, I am adding a new layer to the previous graphic.

g +
geom_point(data = test, aes(x = X, y = Y, color = color)) +
scale_color_identity(labels = "2001 wolves", guide ="legend",
                     name = NULL) -> gg