Maps of single US state illustrating state border, county borders, and terrain background

Question

I wish to generate a map of the state of Kansas in R that displays the border of the state, the borders of each of its counties, and a terrain background. To that end, I use the libraries ggmap and mapdata.

library(ggmap)
library(ggthemes)
library(dplyr)
library(mapdata)

ks_state_map_df <- map_data("state", region="kansas")
ks_counties_map_df <- map_data("county") %>% filter(region=="kansas")
ks_cities_map_df <- us.cities %>% filter(country.etc=="KS")
ks_cities_map_df$name <- gsub(" KS","", ks_cities_map_df$name)
KS_area <- c(left=-102.051744, right=-94.588413, bottom=36.993016, top=40.003162)
KS_map <- get_stadiamap(bbox=KS_area, maptype='outdoors')

ggmap(KS_map) +
  theme_map() +
  geom_polygon(data=ks_state_map_df, mapping=aes(x=long, y=lat, group=group), color="black", fill=NA) + 
  geom_polygon(data=ks_counties_map_df,  mapping=aes(x=long, y=lat), color="gray", fill=NA) + 
  geom_point(data=ks_cities_map_df, mapping=aes(x=long, y=lat)) +
  geom_text(data=ks_cities_map_df, mapping=aes(x=long, y=lat, label=name), hjust=0, nudge_x=0.1)

However, it seems that my above code attempt does not render the polygons correctly (something to do with the group attribute, I believe). How do I need to fix my code?

Edit 1: The terrain map areas should be increased to: KS_area <- c(left=-102.2, right=-94.4, bottom=37.2, top=40.2)

Answer 1

Updated answer

As @r2evans mentions, you need to group each individual feature. But when dealing with geospatial point, line, and polygon data, the sf package is almost always the most robust approach. That's because it allows you set the CRS of each layer etc. I haven't used st_crs() and/or st_transform() in this example and instead am using inherit.aes = FALSE for the sf objects.

Also, I have used the ggrepel package for dealing with overlapping text labels. Getting the geom_text_repel() parameters right is more 'art' than science, but it can be an effective approach.

library(ggmap)
library(ggthemes)
library(dplyr)
library(mapdata)
library(sf)
library(ggrepel)

# sf versions
ks_state_map_df <- map_data("state", region="kansas") %>%
  st_as_sf(coords = c("long","lat")) %>%
  group_by(subregion) %>%
  summarise(geometry = st_combine(geometry)) %>%
  st_cast("POLYGON")

ks_counties_map_df <- map_data("county") %>% 
  filter(region=="kansas") %>%
  st_as_sf(coords = c("long","lat")) %>%
  group_by(subregion) %>%
  summarise(geometry = st_combine(geometry)) %>%
  st_cast("POLYGON")

ks_cities_map_df <- us.cities %>% 
  filter(country.etc=="KS") %>%
  st_as_sf(coords = c("long","lat"))

ks_cities_map_df$name <- gsub(" KS","", ks_cities_map_df$name)

KS_area <- c(left=-102.051744, right=-94.588413, bottom=36.993016, top=40.003162)
KS_map <- get_stadiamap(bbox=KS_area, maptype='outdoors')

ggmap(KS_map) +
  geom_sf(data=ks_state_map_df,
          color="red", # For emphasis
          fill=NA,
          linewidth=2,
          inherit.aes = FALSE) + 
  geom_sf(data=ks_counties_map_df,
          color="gray",
          fill=NA,
          linewidth = 1,
          inherit.aes = FALSE) + 
  geom_sf(data=ks_cities_map_df,
          inherit.aes = FALSE) +
  geom_text_repel(data=ks_cities_map_df,
                  aes(x=st_coordinates(ks_cities_map_df)[,1], 
                      y=st_coordinates(ks_cities_map_df)[,2],
                      label = name),
                  size = 3,
                  max.overlaps = Inf,
                  force = 2) +
  coord_sf(clip = 'off')

Answer 2

You are missing group=group from ks_counties_map_df.

I'll exaggerate the lines for the four geoms just to highlight what is "new".

ggmap(KS_map) +
  theme_map() +
  geom_polygon(data=ks_state_map_df, mapping=aes(x=long, y=lat, group=group), color="black", fill=NA, linewidth=10) + 
  geom_polygon(data=ks_counties_map_df,  mapping=aes(x=long, y=lat, group=group), color="gray", fill=NA, linewidth=10) + 
  geom_point(data=ks_cities_map_df, mapping=aes(x=long, y=lat), size=10) +
  geom_text(data=ks_cities_map_df, mapping=aes(x=long, y=lat, label=name), hjust=0, nudge_x=0.1, size=10)

I'm not entirely sure about ks_state_map_df, it's causing the two black diagonals.