7

I am having clipping problems when I try to combine ggmap with shape files. The example in Kahle and Wickham (2013: 158) works fine because the raster image from ggmap covers the entire shape file. Below is an example of what happens when I try to plot the shape file for U.S. states on a ggmap plot that covers a smaller area. The ggmap shows New York City and I want to overlay it with the borders for U.S. states (just as an example). The resulting map doesn't make any sense. The problem is that the shape file gets clipped and ggplot connects the unclipped points. Below is the code. The shape file is from here. I am just showing the last plot here.

How can I solve this problem?

path <- "PATH TO SHAPEFILE"
library("ggmap")
library("rgdal")

# shapefile
states <- readOGR(dsn = path, layer = "states")
states_df <- fortify(states)
# plot shapefile
plot(states, lwd = 0.1)
ggplot(states_df, aes(long, lat, group = group)) +
    geom_polygon(colour = "black", fill = NA, size = 0.1)


# combine ggmap with shapefile
map <- get_map("new york city", zoom = 10, source = "stamen")
ggmap(map, extent = "device")

ggmap(map, extent = "device") +
    geom_polygon(aes(long, lat, group=group), data = states_df, colour = "red", fill = NA, size = 1)

Kahle, David and Hadley Wickham. 2013. “Ggmap: Spatial Visualization with ggplot2.” The R Journal 5(1):144–61.

enter image description here

user2503795
  • 4,035
  • 2
  • 34
  • 49
  • Is the problem related to this question here: http://stackoverflow.com/questions/13469566/polygons-nicely-cropping-ggplot2-ggmap-at-different-zoom-levels ? – user1965813 May 08 '15 at 12:59

2 Answers2

5

Here is my attempt. I often use GADM shapefiles, which you can directly import using the raster package. I subsetted the shape file for NY, NJ and CT. You may not have to do this in the end, but it is probably better to reduce the amount of data. When I drew the map, ggplot automatically removed data points which stay outside of the bbox of the ggmap image. Therefore, I did not have to do any additional work. I am not sure which shapefile you used. But, GADM's data seem to work well with ggmap images. Hope this helps you.

library(raster)
library(rgdal)
library(rgeos)
library(ggplot2)

### Get data (shapefile)
us <- getData("GADM", country = "US", level = 1)

### Select NY and NJ
states <- subset(us, NAME_1 %in% c("New York", "New Jersey", "Connecticut"))

### SPDF to DF
map <- fortify(states)

## Get a map
mymap <- get_map("new york city", zoom = 10, source = "stamen")


ggmap(mymap) +
geom_map(data = map, map = map, aes(x = long, y = lat, map_id = id, group = group))

enter image description here

If you just want lines, the following would be what you are after.

ggmap(mymap) +
geom_path(data = map, aes(x = long, y = lat, group = group))

enter image description here

jazzurro
  • 23,179
  • 35
  • 66
  • 76
  • Nice results, the one problem you are going to run into using geom_path for polygons, is that if the number of vertexes in the shape file aren't sufficiently high then there is going to be truncation near the edges of the map. – Jonathan Lisic May 08 '15 at 01:37
  • Does this solve the truncation problem `ggmap(mymap) + geom_map(data = map, map = map, aes(x = long, y = lat, map_id = id, group = group), color = "black", fill=NA, size = 0.5)`? I see the problem at the top middle end of the border line but it doesn't seem to be the case with `geom_map` and `fill=NA`. – user2503795 May 08 '15 at 23:13
  • @user2503795 It seems that your code is working well. I am happy for you. :) – jazzurro May 09 '15 at 04:44
  • @jazzurro: Thanks! I think the simplest answer to my question based on your solution is: Use `geom_map` instead of `geom_polygon` and adjust arguments (add `map = [fortified dataframe]` and `map_id = id` as aesthetic to `aes` call). – user2503795 May 10 '15 at 18:46
  • @user2503795 I am very glad to hear that you found your solution based on this answer. Thanks for sharing your own idea here. This will help some other users in the future. :) – jazzurro May 10 '15 at 22:22
  • Thanks! How would you do the opposite? Making what is outside the polygon in transparent? – Agus camacho May 25 '23 at 19:30
2

I would check out this answer, it seems that ggmap as you expected doesn't handle polygon's in an ideal way when you zoom in, namely items not on the plot get truncated causing 'interesting' results with respect to the shape files.

Polygons nicely cropping ggplot2/ggmap at different zoom levels

# transform for good measure
states <- spTransform(states,CRS("+datum=WGS84 +proj=longlat") )

# combine ggmap with shapefile
states_df <- fortify(states)

# get your map
map <-get_map("new york city", zoom = 10, source = "stamen")

a <- ggmap(map, # this is where we get our raster
       base_layer=ggplot(aes(x=long, y=lat), data=states_df), # this defines the region where things are plotted
       extent = "normal",  # this won't work with device, you need normal (see examples in ggmap documentation)
       maprange=FALSE
       ) +
coord_map( # use map's bounding box to setup the 'viewport' we want to see
  projection="mercator",
  xlim= c(attr(map, "bb")$ll.lon, attr(map, "bb")$ur.lon),
  ylim=c(attr(map, "bb")$ll.lat, attr(map, "bb")$ur.lat)
) +
geom_polygon( # plot the polygon
  aes(x=long, y=lat,group=group), data =states_df, color = "red", fill=NA, size = 1)

print(a)

With output: enter image description here

As a side note you might want to check out using the U.S. Census data for state maps, they seem to be of higher quality than the ESRI data set.

ftp://ftp2.census.gov/geo/pvs/tiger2010st/tl_2010_us_state10.zip

As a final note, there are issues with ggmap near the poles so I would also subset your data by the states you are interested in.

Community
  • 1
  • 1
Jonathan Lisic
  • 1,710
  • 1
  • 13
  • 21