facet wrap distorts state maps in R

This question is an extension of a previous question I asked, see Mapping different states in R using facet wrap. I want to display points in the maps, colored by their value. I can get this to work using the following code, which produces a decent looking plot, but I like the look of the 'sepStates' output that Mark Peterson provided better.

mn <- min(test$value); hist(test$value) 
mx <- max(test$value)
ggplot(map_data('state',region=states), aes(x=long,y=lat,group=group)) +
geom_polygon(colour='black',fill='white') + geom_point(data=test,aes
(x=Lon,   y=Lat, group=group, color=value),size=2.5, pch=16) +
facet_wrap(~region, scales = "free", ncol=3) + 
scale_colour_gradientn(limits = c(mn,mx),breaks=c
(2,3,4,5,6,7,8,9),colours=rev(rainbow(7)))

Let's say there's a temperature column added to the zip code data Mark provided.

datatest <- structure(list(zip =c 
("5246", "85118", "85340", "34958", "33022", 
"32716", "49815", "48069", "48551", "58076", "58213", "58524", 
"73185", "74073", "73148", "98668", "98271", "98290"), 
city = c("Chandler","Gold Canyon", "Litchfield
Park", "JensenBeach", "Hollywood", "Altamonte
Springs", "Channing", "Pleasant Ridge", "Flint", "Wahpeton", 
"Ardoch", "Braddock", "Oklahoma City", "Sperry", "Oklahoma City", 
"Vancouver", "Marysville", "Snohomish"), state = c
("AZ", "AZ","AZ", "FL", "FL", "FL", "MI", "MI", "MI", "ND", "ND", "ND",
"OK","OK", "OK", "WA", "WA", "WA"), latitude = c(33.276539, 33.34,33.50835, 
27.242402, 26.013368, 28.744752, 46.186913, 42.472235, 
42.978995, 46.271839, 48.204374, 46.596608, 35.551409, 36.306323, 
35.551409, 45.801586, 48.093129, 47.930902), longitude = c(-112.18717,    
-111.42, -112.40523, -80.224613, -80.144217, -81.22328, -88.04546, 
-83.14051, -83.713124, -96.608142, -97.30774, -100.09497, -97.407537, 
-96.02081, -97.407537, -122.520347, -122.21614, -122.03976),
temperature = c(45,87,33,66,12,69,45,78,23,39,41,104,50,53,40,88,56,29)), 
.Names = c("zip", "city", "state", "latitude", "longitude","temperature"), 
row.names = c(NA, -18L), class = c("tbl_df", "tbl", "data.frame"))

If I want to display this in facets by state with only 1 legend on the right, like the first code above produces, how do I do this? My current code looks like this:

sepStates <- lapply(states, function(thisState){
ggplot(map_data('state',region=thisState), aes(x=long,y=lat,group=group)) +
geom_polygon(colour='black',fill='white') + geom_point(data=datatest
[datatest$region == tolower(thisState),],aes(x=long, y=lat, group=group,
color=temperature),size=2.5, pch=16) + ggtitle(thisState) + coord_map() +
theme_void() + theme(plot.title = element_text(hjust = 0.5))})

plot_grid(plotlist = sepStates)

Using the 'cowplot' package still leads to the states being distorted using the following code. Any help would be greatly appreciated.

legend <- get_legend(ok_plot)
prow <- plot_grid(ok_plot + theme(legend.position="none"),
az_plot + theme(legend.position="none"),
mi_plot + theme(legend.position="right"),
nd_plot + theme(legend.position="none"),
fl_plot + theme(legend.position="none"),
wa_plot + theme(legend.position="right"),
align = 'vh',
hjust = -1,
nrow = 2)

Thank you!

Your approach looks like it is on the right track, but you appear to have forgotten to include a few key components (e.g., the data and where you created the objects ok_plot etc.).

First, I find it easiest to create an object to hold whatever color scale you want to use. I am not sure where you are going with the rainbow palette, as the breaks don't seem to match the extent of your data. So, instead, I am going to use viridis to make one (the palette is great for sequential data).

library(viridis)
sharedScale <-
  scale_color_viridis(limits = range(datatest$temperature))

Then, create your separate plots as before. Note that there is a color mapping added, along with the sharedScale, and that I am suppressing the legends right away.

sepStates <-
  lapply(states, function(thisState){
    ggplot(map_data('state',region=thisState)
           , aes(x=long,y=lat,group=group)) +
      geom_polygon(colour='black',fill='white') +
      geom_point(aes(color = temperature)
                 , data = datatest[datatest$region == tolower(thisState), ]
                 , size = 4
                 , show.legend = FALSE) +
      sharedScale +
      ggtitle(thisState) +
      coord_map() +
      theme_void() +
      theme(plot.title = element_text(hjust = 0.5))
  })

plotting this, as before, will give the maps you want (now with color). However, it won't have a legend. For that, we want to extract a single legend to share across all of the plots. Make sure that the values are the same in every plot (e.g., that you set the limits) when you do this.

sharedLegend <-
  get_legend(
    ggplot(datatest
           , aes(x=long,y=lat,color = temperature)) +
      geom_point() +
      sharedScale +
      theme_void()
  )

Now, we can stitch it all together using nested plot_grid. Inside, we generate the same layout as before (just the six states). Outside, we combine that with the legend. You will have to tinker with rel_widths to match the aspect ratio (and titling) that you want in your final product.

plot_grid(
  plot_grid(plotlist = sepStates)
  , sharedLegend
  , nrow = 1
  , rel_widths = c(6, 1)
)

Gives

Similarly, you may want to play with the legend title/size (or consider putting it under the map instead of next to it).

flatLegend <-
  get_legend(
    ggplot(datatest
           , aes(x=long,y=lat,color = temperature)) +
      geom_point() +
      sharedScale +
      labs(color = "Temp") +
      theme_void() +
      theme(
        legend.title = element_text(size = 18, face = "bold")
        , legend.text = element_text(size = 12)
        , legend.key.size = unit(24, "points")
        , legend.position = "bottom")
  )

plot_grid(
  plot_grid(plotlist = sepStates)
  , flatLegend
  , ncol = 1
  , rel_heights = c(6, 1)
)

Gives