replication.R

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

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_character(),
                 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") %>%
  mutate(origintime = ymd_hms(origintime, tz="UTC")) %>% 
  filter(magnitude > 5) %>% glimpse()
mean_lat <- mean(eq$latitude)
mean_lon <- mean(eq$longitude)

time_seq <- seq.POSIXt(from = ymd_hms("2011-3-1 00:00:00", tz="UTC"),
                       to = ymd_hms("2021-3-1 00:00:00", tz="UTC"),
                       by = "5 mins")
earthdf <- data.frame(date = time_seq,
                      lat = rep(mean_lat, length(time_seq)),
                      lon = rep(mean_lon, length(time_seq)))
sun <- getSunlightPosition(data=earthdf, keep=c("altitude","azimuth"))

# if azimuth < 0 the body is in the east 

sunny = sun %>% 
  select(date, altitude, azimuth) %>%
  mutate(azgroup = ntile(azimuth, 120),
         algroup = ntile(altitude, 90)) %>%
  group_by(azgroup, algroup) %>%
  mutate(minaz = min(azimuth), maxaz = max(azimuth),
         minal = min(altitude), maxal = max(altitude)) %>%
  ungroup()

all <- eq %>% 
  mutate(date = round_date(origintime, unit="5 mins"),
         eq=1, eqI = row_number()) %>%
  select(date, magnitude, eq, eqI) %>% 
  right_join(sunny, by = "date") %>%
  mutate(eq = ifelse(is.na(eq),0,1)) %>%
  group_by(minaz, maxaz, minal, maxal) %>%
  summarise(eq_sum= sum(eq), times = n()) %>%
  ungroup() %>%
  mutate(expected = times * sum(eq_sum)/sum(times),
         poisson = ppois(eq_sum, expected),
         sd = case_when(poisson < 0.95 ~ "normal",
                        TRUE ~ "2+ SD"),
         sd= factor(sd, levels=c("normal", "2+ SD")))
  
  
  
grf <- ggplot(all, aes(xmin=minaz, xmax=maxaz, ymin=minal, ymax=maxal,
                fill=sd, colour=sd)) +
  annotate("rect", xmin=-pi,xmax=pi,ymin=-pi/2,ymax=pi/2,fill="#FFFFFF") +
  annotate("rect", xmin=-pi,xmax=pi,ymin=-pi/2,ymax=0,fill="#DDDDDD") +
  annotate("text", x=-pi, y=.12,label="Day", hjust=0,vjust=0) +
  annotate("text", x=-pi, y=-.12,label="Night", hjust=0,vjust=1) +
  geom_rect() + scale_fill_viridis_d(direction = -1,
                                     name="Frequency\n(poisson\nstandard\ndeviation)\nrelative to\nexpected") + 
  scale_colour_viridis_d(direction = -1,
                         name="Frequency\n(poisson\nstandard\ndeviation)\nrelative to\nexpected") +
  theme_tufte() + 
  geom_hline(yintercept = 0, colour="#000000", size=0.5, linetype=1) +
  geom_hline(yintercept = 0, colour="#FFFFFF", size=0.5, linetype=2) +
  labs(title = "Rate of all 289 New Zealand earthquakes in 10 years with Magnitude > 5
relative to sun direction and height at event time
black/white dashed line represents the horizon",
caption = "Source: Geonet Earthquake Catalogue 2011-03-01 to 2021-03-01",
x="\nDirection of Sun", y="Height of sun (degrees)") +
    scale_x_continuous(breaks = c(-pi+0.006, -pi/2, 0, pi/2, pi),
                       labels = c("N", "E", "S","W","N"))  +
    scale_y_continuous(breaks = c(-pi/4, 0, pi/4),
                       labels = c(-45,0,45)) + 
  theme(plot.background = element_rect(fill="#FFFFFF"),
        panel.background = element_rect(fill = "#FAFAFA", colour = "#FAFAFA"),
        panel.spacing = unit(1.5, "lines"))
grf

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