maps.Rmd

# Stamen maps with ggmap

Here is an example to get started with `ggmap` using `get_stamenmap()` to plot the longitude/latitude maps. The data for the following plots is available at https://simplemaps.com/data/us-cities. The get_stamenmap() function reqiures a bounding box, i.e the top, bottom, left and right latitude/longitude of the map you want to plot. For example, the latitude/longitude for US map are as follows:

```{r bounding_box}
bbox <- c(bottom = 25.75, top = 49 , right = -67, left = -125)
```

You can find these values from https://www.openstreetmap.org. The other important parameters of this function are `zoom` and `maptype`. Higher the zoom level, the more detailed your plot will be. Beaware that ggmap connects to Stamen Map server to download the map, so if your bounding box is large and zoom level is high, it will have to download a lot of data and may take some time. There are differnt types of plots available via Stamen Map like terrain, watercolor, toner which can be set to maptype parameter according to your preference. You can find about avaiable options in help (`?get_stamenmap`). For the following examples the `maptype` is set to ‘toner-lite’.

Let’s plot the US map.

```{r US_map, message=FALSE, warning=FALSE}
library(ggmap)
usmap <- get_stamenmap(bbox = bbox, zoom = 6, maptype = 'toner-lite') 
ggmap(usmap)
```

Great! We have the US map, now let’s use the US population data to see the population density across nation. Notice that we haven’t included Alaska in the map and hence will be removing the data from Alaska.

```{r data, message=FALSE, warning=FALSE}
library(dplyr)
df <- read.csv(unz('resources/ggmap/data/uscities.zip', 'uscities.csv'))
# Removing data of Alaska from dataset
df <- df %>% filter(state_name != 'Alaska')
```

```{r points, message=FALSE, warning=FALSE}
# Population density across US using points
ggmap(usmap) + 
  geom_point(data = df, 
             mapping = aes(x = lng, y = lat, color = population)) +
  ggtitle('Population density across US')
```

This is not good! Most of the points are overlapping and thus it is not easy to interpret what’s going on here. Let’s try alpha blending.

```{r alpha, message=FALSE, warning=FALSE}
# Population density across US using points
ggmap(usmap) + 
  geom_point(data = df, 
             mapping = aes(x = lng, y = lat, color = population),
             stroke= 0, alpha = 0.1) +
  ggtitle('Population density across US')
```

That’s much better! We can now easily identify the areas where population density is more. You might have noticed there is no light blue dot visible on the plot. This is because it must be lying somewhere between those dense areas. One such location is New York, you can find this out by zooming the plot.

***

We can also look at popluation density using `geom_density` as follows

```{r density, message=FALSE, warning=FALSE}
# Population density across US using Density_2d
ggmap(usmap) + 
  geom_density_2d(data = df, 
                  mapping = aes(x = lng, y = lat, color = population)) +
  ggtitle('Population density across US')
```

***

## Mutilayerd plots with ggmaps

We can add multiple layers to the plot as described in earlier chapters. Let’s look at the location of military stations located across US along with population density.

```{r multilayer, message=FALSE, warning=FALSE}
# Location of Military units 
df1 <- df %>% filter(military == TRUE) 

ggmap(usmap) + 
  geom_point(data = df, 
             mapping = aes(x = lng, y = lat, color = population, text = city),
             show.legend = F,
             stroke= 0, alpha = 0.1) +
  geom_point(data = df1, 
             mapping = aes(x = lng, y = lat , text = city), 
             show.legend = F,
             color = 'red') + 
  ggtitle('Military stations across US')
```

***

Let’s zoom the map for state of California.

```{r CA, message=FALSE, warning=FALSE}
# California Boundaries
CAbox <- c(bottom = 32.213, top = 42.163 , right = -113.95, left = -124.585)
camap <- get_stamenmap(bbox = CAbox, zoom = 6, maptype = 'toner-lite')

df3 <- df %>% filter(state_name == 'California')

ggmap(camap) + 
  geom_point(data = df3, 
             mapping = aes(x = lng, y = lat, color = population),
             stroke= 0, alpha = 0.1) + 
  ggtitle('Population density for California')
```

***

## Getting Deeper

This was just a glimpse of what you can do with `ggmaps` using the `get_stamenmap()`. Note that Stamen Maps is not limited to US and can be used to plot any part of the world. If you liked this alternative to Google Maps API, I highly recommend you to check the Stamen Maps website http://maps.stamen.com for more details.