Volcano-City-Connections.Rmd

---
title: "Notebook for connecting volcanos to endangered population centers"
author: "Anne Marshall"
output:
  pdf_document: default
  html_document: default
---

```{r load packages, message = FALSE}

library(ggplot2)
library(patchwork)

install.packages('geosphere')
library(geosphere)
library(dplyr)
library(stringr)
meters_per_mile = 1609.34
```




```{r load the volcano locations}


vol_locations <- read.csv('data/Holocene_volcanos.csv',sep=",",skip = 1)

vol_locations

```

```{r get ranked airports}
vol_locations['Long_rang'] = 1.0/cos(vol_locations['Longitude'])
vol_locations
```

```{r Read population centers}

pop_centers <- read.csv('data/geonames-all-cities-with-a-population-1000-reduced.csv',sep=",") %>%
  rename(geoId = 'ï..geonameId')

pop_locs <- matrix(as.numeric(str_split(pop_centers$Coordinates,",",simplify=TRUE)),ncol=2)
pop_centers["lat"]  = pop_locs[,1]

pop_centers["lon"]  = pop_locs[,2]
pop_centers

```


```{r join population centers with volcanos}

distance_threshold_km = 200
close_cities = data.frame(matrix(ncol=17,nrow=0))

for (irow in 1:nrow(vol_locations)){
  volc <- vol_locations[irow,c("Volcano.Number","Volcano.Name","Country", "Primary.Volcano.Type","Developing.Nation","Longitude","Latitude","Long_rang")]

  print(volc$Volcano.Number)
  # filter candiates for measure to within one deg lat, and close in long (depending on lat)
  candidates <- pop_centers[abs(pop_centers$lat-volc$Latitude)< 1 & abs(pop_centers$lon - volc$Longitude) < abs(as.numeric(volc$Long_rang)),]
  if (nrow(candidates) > 0){
   # calculate geodesic distance
    candidates['pop_vol_dist_km'] <- distm (candidates[,c('lon','lat')],volc[,c('Longitude','Latitude')],fun=distGeo)/1000
    candidates['pop_vol_id'] <- volc['Volcano.Number']

    print(nrow(candidates))
    # filter to cities within range
    match <- candidates[candidates$pop_vol_dist_km < distance_threshold_km,]
  
    print(nrow(match))
    # join back with volcano data
    match <- merge(x=match, y = volc, by.x = 'pop_vol_id',by.y='Volcano.Number')
  
   # append to data file
    close_cities <- rbind(close_cities,match)
  } else {
    print("isolated volcano")
  }
  print(nrow(close_cities))
  }

#Save output

write.csv(x=close_cities,file="data/volcano_cities.csv",row.names=FALSE)

```

```{r plot world}

ggplot() +
  geom_map(
    data = world, map = world,
    aes(long, lat, map_id = region)
  ) 
```



```{r hacky stuff to make the NYC plot}


#for(airport_num in 1:nrow(ranked_airports_exp)){

counties <- map_data("county")

#%>%
#  filter(group == 166 | group ==176 | group ==172 | group ==178 | group ==180 | group ==210) 

chosen <- counties %>%
  filter(group == 1825)

 candidate_counties
 
 plot_counties  <- map_data("county")%>%
  filter(group == 1835 | group ==1818 | group == 1824 | group == 1825 | group == 1797 | group == 1750) 


plain_theme = theme(axis.text=element_blank()) + theme(panel.background = element_blank(), panel.grid = element_blank(), axis.ticks = element_blank())


ap2 <- ap %>%
  rename(lat=Latitude) %>%
  rename(long=Longitude)

# Add Laguardia for comparison
	LGA	
40.7772
-73.8726

lga <-data.frame(IATA_Code="LGA", lat=40.7772,long=-73.8726)




 plot_counties_center <- plot_counties %>%
   group_by(group) %>%
   summarize(lat=mean(lat),long=mean(long)) %>%
  mutate(name=c("Hudson, NJ","Bronx","Kings","Nassau","New York","Queens")) %>%
    mutate(casrn= c(34017L,36005L,36047L,36059L, 36061L,36081L)) 
 

 plot_counties_center
 
 pcounties <- merge(plot_counties_center,candidate_counties,by="casrn") %>%
   mutate(pop_text = round(DaytimePopDensity/100)*1000)

 pcounties
 
 #**I am trying to populations into millions right now!***
 
 png("NYC.png",width=6,height=4,units="in",res=1200)
 ggplot() + 
  aes(x=long, y = lat) +
  geom_polygon(data= plot_counties, aes(group=group), fill="grey90",color="grey",size=.3) +
  geom_polygon(data= chosen, aes(group=group), fill="skyblue2",color="grey",size=.3) +
  geom_text(data=ap2,aes(label=IATA_Code),nudge_y=.01)+
  geom_point(data=ap2,  color = "black",size=1) +
  geom_point(data=ap2,  color = "red",size=.5) +
   geom_text(data=lga,aes(label=IATA_Code),nudge_y=.01)+
  geom_point(data=lga,  color = "black",size=1) +
  geom_point(data=lga,  color = "blue",size=.5) +
  geom_text(data=pcounties,aes(label=name.x),nudge_y=.022,size=4) +
   geom_text(data=pcounties,aes(label=pop_text),nudge_y=.005,size=3) +
  labs(x="", y="") +
  coord_fixed(1.3) +
  plain_theme

 dev.off()

  
```