maps-ggplot

library(tidyverse)

── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
✔ ggplot2 3.3.6      ✔ purrr   0.3.4 
✔ tibble  3.1.8      ✔ dplyr   1.0.10
✔ tidyr   1.2.1      ✔ stringr 1.4.1 
✔ readr   2.1.2      ✔ forcats 0.5.2 
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

library(mdsr)
library(sf)

Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE

library(ggspatial)
library(leaflet)
library(tidygeocoder)
library(mapproj)

Loading required package: maps

Attaching package: 'maps'

The following object is masked from 'package:purrr':

    map

library(maps)

CholeraDeaths_Book <-mdsr::CholeraDeaths
st_crs(CholeraDeaths_Book) = st_crs(CholeraDeaths_Book) 

old-style crs object detected; please recreate object with a recent sf::st_crs()
old-style crs object detected; please recreate object with a recent sf::st_crs()

Let’s find the latitude and longitude for a couple of places.

addresses<- tibble(address=c("50 Payson Ave, Easthampton, MA","Holyoke Community College, Holyoke,MA")) |>
geocode(address, method="osm") |>
  mutate(name = c("City Hall","HCC"))

Passing 2 addresses to the Nominatim single address geocoder

Query completed in: 2 seconds

addresses

# A tibble: 2 × 4
  address                                 lat  long name     
  <chr>                                 <dbl> <dbl> <chr>    
1 50 Payson Ave, Easthampton, MA         42.3 -72.7 City Hall
2 Holyoke Community College, Holyoke,MA  42.2 -72.7 HCC      

Let’s plot those with ggplot.

ggplot(addresses) +
  geom_point(aes(long,lat)) 

  #coord_quickmap()

Not very exciting. We’ll add it it in a moment.

Let’s build a map of Massachusetts with the map_data() function

ma_counties <- map_data("county", "massachusetts") %>% 
  select( long, lat, group, id = subregion)

head(ma_counties)

       long      lat group         id
1 -70.67435 41.73997     1 barnstable
2 -70.53683 41.79727     1 barnstable
3 -70.53683 41.79727     1 barnstable
4 -70.51392 41.78008     1 barnstable
5 -70.47954 41.75716     1 barnstable
6 -70.41078 41.73425     1 barnstable

Notice there are a lot of points for barnstable. Those are the vertices of a polygon. Let’s draw the polygon with the points.

ma_counties |> ggplot(aes(long,lat,group=group))+
  geom_polygon(fill = "white", colour = "grey50")+
  geom_point()+
  coord_quickmap()

Let’s grab just Massachusetts cities from the data frame us.cities.

ma_cities <- us.cities %>% 
  filter(country.etc =="MA") |>
  select(name, long, lat)

head(ma_cities)

                name   long   lat
1         Andover MA -71.14 42.65
2       Arlington MA -71.16 42.42
3       Attleboro MA -71.30 41.93
4 Barnstable Town MA -70.30 41.70
5         Beverly MA -70.84 42.56
6       Billerica MA -71.26 42.56

Exercise: Plot MA cities onto our county map.

ma_counties |> ggplot(aes(long,lat))+
  geom_polygon(fill = "white", colour = "grey50",aes(group=group)) +
  #geom_point(data=ma_cities,.......)
  
  coord_quickmap()

class(ma_counties)

[1] "data.frame"

class(ma_cities)

[1] "data.frame"

Let’s get some more accurate government data.

Go here and download some shape files.

Saving them will be a bit of work. First we need to find those files.

states <- read_sf("~/Documents/MTH_190_Intro_to_DS/course-materials/cb_2021_us_state_500k") |>
  sf::st_transform('+proj=longlat +datum=WGS84')

counties <- read_sf("~/Documents/MTH_190_Intro_to_DS/course-materials/cb_2021_us_county_5m") |>
  st_transform('+proj=longlat +datum=WGS84')

leaflet(states) |>
  addTiles() |>
  addPolygons()

Just_MA <-states |> 
  filter(STUSPS =="MA") 
  
leaflet(Just_MA) |>
  addTiles() |>
  addPolygons()

Just_MA_counties <- counties |> filter(STUSPS == "MA")

leaflet() |>
  addTiles() |>
  addPolygons(data = Just_MA_counties) |>
  addMarkers(data = addresses, lng=~long, lat=~lat, popup = ~name)

Just_MA_counties <- counties |> filter(STUSPS == "MA")

leaflet() |>
  addTiles() |>
  addPolygons(data = Just_MA_counties) |>
  addMarkers(data = addresses, lng=~long, lat=~lat, popup = ~name, ) |>
  addCircleMarkers(data = ma_cities, lng=~long, lat=~lat, popup = ~name, color='black',clusterOptions = markerClusterOptions() ) 

leaflet() |>
  addTiles() |>
  addPolygons(data = Just_MA_counties) |>
  addMarkers(data = addresses, lng=~long, lat=~lat, popup = ~name, ) |>
  addCircleMarkers(data = ma_cities, lng=~long, lat=~lat, popup = ~name, color='black',clusterOptions = markerClusterOptions() ) 

Download some more shape files from here.

Exercise: Plot the Places we should care about from MA on our map.

Exercise: Plot the American Indian Tribal Subdivisions.