public
Last active

Shint App mapping earthquakes by magnitude, country and year

  • Download Gist
global.R
R
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
library(maps)
library(shiny)
library(stringr)
 
worlddf <- map("world")
 
# restrict selection to coutries only
locations <- worlddf$names
 
countries <-subset(locations,is.na(str_locate(locations,":")))
countries <- countries[!is.na(countries)]
countryChoice <- sort(countries)
countryChoice <- append("World",countryChoice)
 
# Note this leaves some countries e.g Japan missing
# and includes others e.g. USSR that no longer exist
 
### Potential Enhancements ###
 
# Extend to all current countries
# Look at alternative mapping packages
# Extend to shorter time periods as alternatives to years
# Add source giving lower magnitudes to show at smaller geographical areas e.g. counties
# Assess other plotting packages e.g. ggplot2
# Offer user option of splitting magnitudes by count rather than equi-.... cuts
# Add more information on graph and/or table on time position of some/all quakes
server.R
R
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84
# Install as Required
 
library(maptools)
library(rgdal)
library(XML)
 
 
 
shinyServer(function(input, output) {
 
# create and plot mapmap
output$map <- reactivePlot(function() {
# collect variables for insertion into url
lmag <- input$mag[1]
umag <- input$mag[2]
syear <- input$year[1]
eyear <- input$year[2]
# obtain and process data into a dataframe
myURL <- paste0("http://neic.usgs.gov/cgi-bin/epic/epic.cgi?SEARCHMETHOD=1&FILEFORMAT=6&SEARCHRANGE=HH&SYEAR=",syear,"&SMONTH=01&SDAY=01&EYEAR=",eyear,"&EMONTH=1&EDAY=31&LMAG=",lmag,"&UMAG=",umag,"&NDEP1=&NDEP2=&IO1=&IO2=&CLAT=0.0&CLON=0.0&CRAD=0.0&SUBMIT=Submit+Search")
basicInfo <- htmlParse(myURL, isURL = TRUE)
data <- xpathSApply(basicInfo, "//*/pre/text()", xmlValue)
eq <- read.table(text = data, fill=TRUE, sep = ',',header=TRUE)
# create plot (based on http://statistical-research.com/earthquakes-over-the-past-7-days/
plot.new()
if (input$country=="World") {
my.map <- map("world", interior = FALSE, plot=T)
map("world", boundary = FALSE, col="gray", add = TRUE) # this works after having defined map above
} else {
my.map <- map("world",regions=input$country, interior = FALSE, plot=T)
}
# split magnitude range into four groups
cut(eq$Magnitude,4)
labs <- levels(cut(eq$Magnitude, 4))
ranges <-as.data.frame(cbind(lower = as.numeric( sub("\\((.+),.*", "\\1", labs) ),
upper = as.numeric( sub("[^,]*,([^]]*)\\]", "\\1", labs) )))
# make point size and colour dependent on magnitude level
eq$mag.size <- NULL
eq$mag.size[eq$Magnitude>=ranges$lower[1] & eq$Magnitude<ranges$upper[1]] <- .75
eq$mag.size[eq$Magnitude>=ranges$lower[2] & eq$Magnitude<ranges$upper[2]] <- 1.0
eq$mag.size[eq$Magnitude>=ranges$lower[3] & eq$Magnitude<ranges$upper[3]] <- 1.5
eq$mag.size[eq$Magnitude>=ranges$lower[4]] <- 2.0
eq$mag.col <- NULL
eq$mag.col[eq$Magnitude>=ranges$lower[1] & eq$Magnitude<ranges$upper[1]] <- 'blue'
eq$mag.col[eq$Magnitude>=ranges$lower[2] & eq$Magnitude<ranges$upper[2]] <- 'green'
eq$mag.col[eq$Magnitude>=ranges$lower[3] & eq$Magnitude<ranges$upper[3]] <- 'orange'
eq$mag.col[eq$Magnitude>=ranges$lower[4]] <- 'red'
points(x=eq$Lon,y=eq$Lat,pch=16,cex=eq$mag.size, col=eq$mag.col)
legend1 <- paste0("M ",ranges$lower[1],"-",ranges$upper[1])
legend2 <- paste0("M ",ranges$lower[2],"-",ranges$upper[2])
legend3 <- paste0("M ",ranges$lower[3],"-",ranges$upper[3])
legend4 <- paste0("M ",ranges$lower[4],"+")
legend('bottomright',c(legend1,legend2,legend3,legend4), ncol=2,
pch=16,cex=0.75, col=c('blue','green','orange','red'))
box()
})
# Create Heading
output$caption <- reactiveText(function() {
if (input$year[1]==input$year[2]) {
paste0("Earthquakes between Magnitude ",input$mag[1]," and ",input$mag[2],", ",input$year[1])
} else {
paste0("Earthquakes between Magnitude ",input$mag[1]," and ",input$mag[2],", ",input$year[1]," - ",input$year[2])
}
})
})
ui.R
R
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51
shinyUI(pageWithSidebar(
# Application title
headerPanel("Earthquake Maps POC"),
# Sidebar with information, controls to select time period, magnitude and geographical coverage
sidebarPanel(
helpText(
p("A Proof-of-Concept app requiring more work as shown in global.r file of",
a("associated gist.", href="https://gist.github.com/4208051"),
 
"Inspired by",
a("this blog post", href="http://statistical-research.com/earthquakes-over-the-past-7-days/?utm_source=rss&utm_medium=rss&utm_campaign=earthquakes-over-the-past-7-days")
) ,
p(
"Select a time period, magnitude and region to map. For guidance in plotting, there
are around 1.5 million annually worldwide with a magnitude of 2 or greater but only 150, 6+"),
p("According to the NEIC. the catalog contains data for M2.5 and greater U.S. earthquakes and M4.5 and greater for Worldwide earthquakes
although some smaller values do appear")
),
wellPanel(
sliderInput(inputId = "mag",
label="Select a range of Magnitudes",
min = 1.7, max = 10, step = 0.1, value = c(4,10)),
# radioButtons(inputId="region", label="Select region US has ",choices=c("World (Mag 4.5+)","USA (Mag 2+)"))
selectInput("country", "Select Country" ,countryChoice ,selected="World"),
sliderInput(inputId = "year",
label="Select a range of Years",
min = 1990, max = 2012, step = 1,
format="####",
ticks=FALSE,
value = c(2011,2012))
)
),
# Show the plot
mainPanel(
h4(textOutput("caption")),
plotOutput("map")
)
))

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.