chart_challenge_2021_day_6.R

library(readr)
library(suncalc)
library(dplyr)
library(ggplot2)
library(ggthemes)
library(lubridate)

quake10yr <- c("https://quakesearch.geonet.org.nz/csv?bbox=165.5420,-49.1817,-178.9893,-32.2871&minmag=0&mindepth=0&startdate=2020-02-01&enddate=2021-03-01T0:00:00",
               "https://quakesearch.geonet.org.nz/csv?bbox=165.5420,-49.1817,-178.9893,-32.2871&minmag=0&mindepth=0&startdate=2019-02-01&enddate=2020-02-01",
               "https://quakesearch.geonet.org.nz/csv?bbox=165.5420,-49.1817,-178.9893,-32.2871&minmag=0&mindepth=0&startdate=2018-01-01&enddate=2019-02-01",
               "https://quakesearch.geonet.org.nz/csv?bbox=165.5420,-49.1817,-178.9893,-32.2871&minmag=0&mindepth=0&startdate=2017-03-01&enddate=2018-01-01",
               "https://quakesearch.geonet.org.nz/csv?bbox=165.5420,-49.1817,-178.9893,-32.2871&minmag=0&mindepth=0&startdate=2016-11-01&enddate=2017-03-01",
               "https://quakesearch.geonet.org.nz/csv?bbox=165.5420,-49.1817,-178.9893,-32.2871&minmag=0&mindepth=0&startdate=2015-11-01&enddate=2016-11-01",
               "https://quakesearch.geonet.org.nz/csv?bbox=165.5420,-49.1817,-178.9893,-32.2871&minmag=0&mindepth=0&startdate=2014-10-01&enddate=2015-11-01",
               "https://quakesearch.geonet.org.nz/csv?bbox=165.5420,-49.1817,-178.9893,-32.2871&minmag=0&mindepth=0&startdate=2013-10-01&enddate=2014-10-01",
               "https://quakesearch.geonet.org.nz/csv?bbox=165.5420,-49.1817,-178.9893,-32.2871&minmag=0&mindepth=0&startdate=2012-11-01&enddate=2013-10-01",
               "https://quakesearch.geonet.org.nz/csv?bbox=165.5420,-49.1817,-178.9893,-32.2871&minmag=0&mindepth=0&startdate=2011-10-01&enddate=2012-11-01",
               "https://quakesearch.geonet.org.nz/csv?bbox=165.5420,-49.1817,-178.9893,-32.2871&minmag=0&mindepth=0&startdate=2011-03-01T0:00:00&enddate=2011-10-01")

if(!file.exists("eq10yrs.csv")){
  q_list <- lapply(quake10yr, read_csv, col_types = cols(.default = col_character()))
  q_df <- bind_rows(q_list)
  write.csv(q_df, file = "eq10yrs.csv",row.names = FALSE)
}

# catalogue in UTC, so is suncalc so it is all good

eq <- read_csv("eq10yrs.csv",
               col_types=cols(
                 .default = col_double(),
                 publicid = col_character(),
                 eventtype = col_character(),
                 origintime = col_datetime(format = ""),
                 modificationtime = col_datetime(format = ""),
                 magnitudetype = col_character(),
                 depthtype = col_character(),
                 evaluationmethod = col_character(),
                 evaluationstatus = col_character(),
                 evaluationmode = col_character(),
                 earthmodel = col_character())) %>%
  filter(is.na(eventtype) | eventtype == "earthquake")
mean_lat <- mean(eq$latitude)
mean_lon <- mean(eq$longitude)
time_seq <- seq.POSIXt(from = ymd_hms("2011-3-1 00:00:00"),
                       to = ymd_hms("2021-3-1 00:00:00"),
                       by = "5 mins")
earthdf <- data.frame(date = time_seq,
                      lat = rep(mean_lat, length(time_seq)),
                      lon = rep(mean_lon, length(time_seq)))
moon <- getMoonPosition(data=earthdf, keep=c("altitude","azimuth"))
sun <- getSunlightPosition(data=earthdf, keep=c("altitude","azimuth"))

# if azimuth < 0 the body is in the east 
# so using that to represent body position as 0 dawn 180 dusk
sunny = sun %>% mutate(ew = ifelse(azimuth < 0, "east", "west"),
  daymins = hour(date) * 60 + minute(date),
  degrees = altitude * 180/pi,
  sun_angle = case_when(ew == "east" & degrees >= 0 ~ degrees,
                        ew == "east" & degrees < 0 ~ 360 - abs(degrees),
                        ew == "west" & degrees > 0 ~ 180 - degrees,
                        ew == "west" & degrees <= 0 ~ 180 + abs(degrees)),
  sun_angle=round(sun_angle)) %>%
  arrange(date)
moony = moon %>% mutate(ew = ifelse(azimuth < 0, "east", "west"),
                       daymins = hour(date) * 60 + minute(date),
                       degrees = altitude * 180/pi,
                       moon_angle = case_when(ew == "east" & degrees >= 0 ~ degrees,
                                             ew == "east" & degrees < 0 ~ 360 - abs(degrees),
                                             ew == "west" & degrees > 0 ~ 180 - degrees,
                                             ew == "west" & degrees <= 0 ~ 180 + abs(degrees)),
                       moon_angle=round(moon_angle)) %>%
  arrange(date)
space <- sunny %>% select(date, sun_angle)
space$moon_angle <- moony$moon_angle

eqs <- eq %>% mutate(date = round_date(origintime, "5 mins")) %>%
  select(date) %>% 
  inner_join(space, by="date") %>%
  count(sun_angle, moon_angle) %>%
  rename(quakes=n)

all <- space %>% 
  count(sun_angle, moon_angle) %>%
  rename(times=n) %>%
  left_join(eqs, by=c("sun_angle","moon_angle")) %>%
  mutate(quakes = ifelse(is.na(quakes), 0, quakes),
         expected = round(times * sum(quakes) / sum(times)),
         rate= case_when(quakes > expected ~ "above average",
                         quakes == expected ~ "average",
                         quakes < expected ~ "below average"),
         rate = factor(rate, levels=c("above average", "average", "below average")))

grf <- all %>% 
  ggplot(aes(x=sun_angle, y=moon_angle, fill=rate)) + coord_fixed() +
  scale_fill_viridis_d(name="Rate") + geom_tile() + theme_tufte() +
  labs(title = "Rate of 223115 New Zealand earthquakes relative\nto daily sun and moon positon at event time.
Angle relative to epicentre 0: east horizon, 90: vertical", x="Sun Angle", y="Moon Angle",
caption = "Source: Geonet Earthquake Catalogue 2011-03-01 to 2021-03-01") + 
  scale_x_continuous(breaks = c(0,90,180,270,360)) +
  scale_y_continuous(breaks = c(0,90,180,270,360)) +
  annotate("segment", x=370, xend=370, y=180, yend=360) +
  annotate("segment", y=370, yend=370, x=180, xend=360) 

ggsave(filename="~/Desktop/ouput.png", plot=grf, dpi=72, units="in",
       height = 5.556, width=9.877)