Shint App mapping earthquakes by magnitude, country and year

global.R

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

# 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

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")
    
    )
  
  
  ))