How to preserve scale_fill_color from first plot using shapefiles in ggplot2

Question

I need help with scale_fill_manual using a shapefile in ggplot2. I have tried many thing and I am finally posting so hopefully someone will be able to give me a hint. I am basically ploting a shapefile and use scale_fill_manual to vizualise it with custom colors, I then overlay some point on it but when I tried to include my new points in the legend, my original colors are there but the values get all messed up. The upper part ploting the shapefile works fine, but the bottom part when overlaying the new points is where I need help. I have some comments inline. See below: The path to download the shapefile is: 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

I get the following graph(good):

All of the above code works as expected. The problem that I am having is down below. The code below works well overlaying the points from a subset of the same shapefile above. However, I lose my scale_fill_manual colors when I try to add the new points to the legend.

#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')

If I try to include '2001 wolves' in my legend, my colors stay in the correct order. However, my legend values get all messed up. I tried to re-arrange them with a different palette but it only makes it worse because the colors and labels don't line up with the corresponding color. I also would like help with removing the dots from the legend. How can I get my colors back to my original Palette1 used above on the original plot? Probably a simple thing but I have spent many hours trying and can't figure it out. Thanks beforehand.

I get this plot. Notice the values are all over. I need the values to be in the same order as on first plot.

EDIT: This is what my plots show behind the scenes. The first plot has the following color order:

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

And my second plot has this color order which is different from the first one. I wonder how I can make the second match the first color order.

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

Have you seen this post? https://stackoverflow.com/questions/11774262/how-to-extract-the-fill-colours-from-a-ggplot-object — mrhellmann, Jan 04 '20 at 03:28
I just took a look at it. This is what I found: Their color order is different:> 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 > I wonder how to make the second color order match the first one. — Salvador, Jan 04 '20 at 04:17
I have been carefully following your code and rewriting it in my own way now. In the final part, you have `2001 wolves`. Where does this come from? I cannot find it in any data you created beforehand. Can you explain it? — jazzurro, Jan 04 '20 at 05:18
Jazzurro: 2001 wolves does no exist anywhere in the code. I created it inside the aes() call to create an entry on the legend manually. If you look at the ggplot documentation it is how one adds a field to an existent legend. — Salvador, Jan 04 '20 at 05:26

jazzurro · Accepted Answer · 2020-01-04T14:04:48.637

Here is what I tried for you. I went through all of your code and had the impression that you made data processing complicated in my view. I used to use the sp approach and write codes like yours. I think this approach made you "twist" your data processing somewhere (e.g. the moment you used merge(), for example). Here I wrote your code in another way in order to deliver the expected outcome. I left explanation in the script below. It seems to me that the takeaway is to avoid some tricky data manipulation. One way of doing it is to use the sf and tidyverse packages. I hope this will help you.

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

Jazzurro; Thanks for showing a different way to achieve the same thing. However, how can you add the blue dots to the legend? — Salvador, Jan 04 '20 at 05:41
@Salvador I am occupied right now. Let me get back to you later. — jazzurro, Jan 04 '20 at 05:58
@ Jazzurro, No problem. Thanks for checking into it. The entry in the legend and re-organizing the legend entries in the correct order is what is keeping me awake. — Salvador, Jan 04 '20 at 06:07
@Salvador I revised the final graphic. I am traveling right now. So I won't be able to help you for some time. I hope this is enough for you to keep going. — jazzurro, Jan 04 '20 at 07:24
This is great, Thank you and have a great trip. I will check it as solved by you. — Salvador, Jan 04 '20 at 07:41
@Salvador I hope you can learn something from my suggestion and make a progress. — jazzurro, Jan 04 '20 at 08:39

score 0 · Answer 2 · answered Jan 04 '20 at 04:01

0

I'm not cool enough to comment, but I'm not confident you need a fill aesthetic on the dots. Adding new values to the fill scale is likely shifting the colors. If deleting it doesn't work, try giving the dots fill=NA inside geom_point() but not in aes().

answered Jan 04 '20 at 04:01

user3312530

1

I tried fill=NA outside the aes() and it gives me the original colors but it deletes the entry from the legend. In order to create a field in the legend for the new points I think it needs to be inside the aes(), however, that changes my color order. – Salvador Jan 04 '20 at 04:42

How to preserve scale_fill_color from first plot using shapefiles in ggplot2

2 Answers2