DAT-4313 - GIS with R

Meteorites and Population Overlay

Author

Dr. V

library(ggplot2)
library(ggthemes)
library(socviz)
library(maps)
library(mapproj)
library(viridis)

GET ARKANSAS POPULATION DATA BY COUNTY DIRECTLY FROM THE US CENSUS VIA API

head(as_tibble(ar_pop))

# A tibble: 6 × 6
  GEOID NAME                   variable estimate   moe                  geometry
  <chr> <chr>                  <chr>       <dbl> <dbl>        <MULTIPOLYGON [°]>
1 05081 Little River County, … B01003_…    12024    NA (((-94.48558 33.65331, -…
2 05121 Randolph County, Arka… B01003_…    18619    NA (((-91.40687 36.49712, -…
3 05013 Calhoun County, Arkan… B01003_…     4773    NA (((-92.77672 33.53926, -…
4 05061 Howard County, Arkans… B01003_…    12779    NA (((-94.2549 34.3462, -94…
5 05099 Nevada County, Arkans… B01003_…     8292    NA (((-93.48322 33.47617, -…
6 05103 Ouachita County, Arka… B01003_…    22606    NA (((-93.11638 33.45245, -…

CREATE A CHOLOPLETH USING LEAFLET FOR ARKANSAS COUNTIES

MapPalette <- colorQuantile(palette = "viridis", domain = ar_pop$estimate, n = 20)

library(sf)
ar_pop %>% 
  st_transform(crs = "+proj=longlat +datum=WGS84") %>% 
  leaflet(width = "100%", height = 500) %>% 
  addProviderTiles(provider = "Esri.WorldPhysical") %>% 
  addPolygons(popup = ~NAME,
              stroke = FALSE,
              smoothFactor = 0,
              fillOpacity = 0.7,
              color = ~ MapPalette(estimate)) %>% 
  addLegend("bottomright", 
            pal = MapPalette,
            values = ~ estimate,
            title = "Population Percentiles",
            opacity = 1)

head(read.csv("Meteorite_Landings.csv"))

      name  id nametype    recclass mass..g. fall year    reclat    reclong
1   Aachen   1    Valid          L5       21 Fell 1880  50.77500    6.08333
2   Aarhus   2    Valid          H6      720 Fell 1951  56.18333   10.23333
3     Abee   6    Valid         EH4   107000 Fell 1952  54.21667 -113.00000
4 Acapulco  10    Valid Acapulcoite     1914 Fell 1976  16.88333  -99.90000
5  Achiras 370    Valid          L6      780 Fell 1902 -33.16667  -64.95000
6 Adhi Kot 379    Valid         EH4     4239 Fell 1919  32.10000   71.80000
        lat       long
1  50.77500    6.08333
2  56.18333   10.23333
3  54.21667 -113.00000
4  16.88333  -99.90000
5 -33.16667  -64.95000
6  32.10000   71.80000

PLOT Meteorite LOCATIONS IN ARKANSAS

# load all USA Starbucks locations
# load("starbucksUS.RData")
# head(starbucksUS)

# filter for just AR
# sbar <- starbucksUS %>% filter(state=="AR")
sbar <- read.csv("Meteorite_Landings.csv") %>% filter(lat >= 33, lat <= 36.5, long >= -94.5, long <= -90)

# dataset is old -- let's add Siloam Springs Starbucks manually:)
# sbar <- sbar %>% 
  # add_row(name="Highway 412", city="Siloam Springs", state="AR", country="US",
          # latitude=36.181498, longitude=-94.507065)

# plot leaflet
sbar %>% leaflet(width = "100%") %>% 
             addTiles() %>% 
             setView(-94.0, 36.0, zoom = 8) %>% 
             addMarkers(lat = ~lat, 
                                 lng = ~long,
                                 icon = makeIcon("asteroid-icon.png", iconWidth = 30, iconHeight = 30))

COMBINE CHOROPLETH AND LOCATIONS – OVERLAY Meteorite LOCATIONS ONTO CHOROPLETH

MapPalette <- colorQuantile(palette = "viridis", domain = ar_pop$estimate, n = 20)


ar_pop %>% 
  st_transform(crs = "+proj=longlat +datum=WGS84") %>% 
  leaflet(width = "100%", height = 500) %>% 
  addProviderTiles(provider = "Esri.WorldPhysical") %>% 
  addPolygons(popup = ~NAME,
              stroke = FALSE,
              smoothFactor = 0,
              fillOpacity = 0.7,
              color = ~ MapPalette(estimate)) %>% 
  addLegend("bottomright", 
            pal = MapPalette,
            values = ~ estimate,
            title = "Population Percentiles",
            opacity = 1) %>% 
addMarkers(data = sbar, 
              lat = sbar$lat,
              lng = sbar$long,
              icon = makeIcon("asteroid-icon.png", iconWidth = 30, iconHeight = 30))