25 Day 25 (April 23)

25.1 Announcements

  • Please send email to and to request 20 min time for presentations between Tuesday April 30 - Thursday May 9.

  • TEVALS will be available soon

    • Your feedback is valued
    • I’ve got lots of questions (e.g., activities and. assignments)

25.2 Earthquake data

  • New York Times article

  • New York Times articles

  • Earthquake data from Kansas

    library(sf)
library(sp)
library(raster)
library(lubridate)
library(animation)    
library(gifski)


# Download shapefile of Kansas from census.gov
download.file("http://www2.census.gov/geo/tiger/GENZ2015/shp/cb_2015_us_state_20m.zip", destfile = "states.zip")
unzip("states.zip")
sf.us <- st_read("cb_2015_us_state_20m.shp",quiet = TRUE)
sf.kansas <- sf.us[48,6]
sf.kansas <- as(sf.kansas, 'Spatial')

# Load earthquake data
url <- "https://www.dropbox.com/scl/fi/uhicc7qo4zcxeq79y1eh9/ks_earthquake_data.csv?rlkey=lbq8kxzx9g067domp46ffh7jw&dl=1"
df.eq <- read.csv(url)
df.eq$Date.time <- ymd_hms(df.eq$Date.time)
pts.eq <- SpatialPointsDataFrame(coords = df.eq[,c(3,2)], data = df.eq[,c(1,4)], proj4string = crs(sf.kansas))

# Plot spatial map earthquake data
par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
plot(sf.kansas,main="")
points(pts.eq,col=rgb(0.4,0.8,0.5,0.3),pch=21,cex=pts.eq$Magnitude/3)
legend("right",inset=c(-0.25,0),legend = c(1,2,3,4,5), bty = "n", text.col = "black", 
   pch=21, cex=1.3,pt.cex=c(1,2,3,4,5)/3,col=rgb(0.4,0.8,0.5,0.6))

    # Plot timeseries of earthquake data
plot(as.numeric(names(table(year(df.eq$Date.time)))),c(table(year(df.eq$Date.time))),
 xlab="Year",ylab="Number of Earthquakes in KS",pch=20)

    plot(df.eq$Date.time,df.eq$Magnitude,pch=20,cex=.2,xlab="Year",ylab="Earthquake magnitude")

    • Animation of Kansas Earthquake data
    # Make animation of earthquake data
date <- seq(as.Date("1977-1-1"), as.Date("2020-10-31"), by = "year")
t <- round(time_length(interval(min(date),pts.eq$Date.time),"year"))

for(i in 1:length(date)){
  par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
  plot(sf.kansas,main=year(date[i]))
  legend("right",inset=c(-0.25,0),legend = c(1,2,3,4,5), bty = "n", text.col = "black", 
     pch=20, cex=1.3,pt.cex=c(1,2,3,4,5)/3,col=rgb(0.4,0.8,0.5,1))
  if(length(which(t==i))>0){
  points(pts.eq[which(t==i),],col=rgb(0.4,0.8,0.5,1),pch=20,cex=pts.eq[which(t==i),]$Magnitude/3)
                        }
                    }

25.3 Spatio-temporal models for earthquake data