library(sp)
library(plotrix)
library(raster)
library(rgeos)
setwd("C:/Visualize Earthquake")
#http://earthquake.usgs.gov/earthquakes/feed/v1.0/csv.php
#From this address we can download csv files with locations of Earthquakes
#For this experiment we will download the summary of the last 30 days
#This dataset is updated every 15 minutes so your output may differ from the following
URL <- "http://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/all_month.csv"
Earthquake_30Days <- read.table(URL, sep = ",", header = T)
#Check the data.frame dimensions
nrow(Earthquake_30Days)
ncol(Earthquake_30Days)
#We can check the contents of the table with the following function
str(Earthquake_30Days)
#The function names() is used to generate a list of the column names of a data.frame
names(Earthquake_30Days)
#Now we can transform the data.frame into a Spatial object
coordinates(Earthquake_30Days)=~longitude+latitude
str(Earthquake_30Days)
#Slot data
str(Earthquake_30Days@data)
#Check and set the projection
Earthquake_30Days@proj4string
#We can set the projection to WGS84 using the following line
projection(Earthquake_30Days)=CRS("+init=epsg:4326")
#http://spatialreference.org/ref/epsg/wgs-84/
#Downloading polygons with the border of each country
#Download, unzip and load the polygon shapefile with the countries' borders
download.file("http://thematicmapping.org/downloads/TM_WORLD_BORDERS_SIMPL-0.3.zip",destfile="TM_WORLD_BORDERS_SIMPL-0.3.zip")
unzip("TM_WORLD_BORDERS_SIMPL-0.3.zip",exdir=getwd())
polygons <- shapefile("TM_WORLD_BORDERS_SIMPL-0.3.shp")
#Let us check the time variable
Earthquake_30Days$time[1]
#The time variable has the following format: year-month-dayThour:minutes:second.milliseconds
#Now that we understood the format we can change this variable from a factor to a time stamp
conv.time <- function(vector){
split1 <- strsplit(paste(vector),"T")
split2 <- strsplit(split1[[1]][2],"Z")
fin <- paste0(split1[[1]][1],split2[[1]][1])
paste(as.POSIXlt(fin,formate="%Y-%m-%d%H:%M:%OS3"))
}
conv.time(Earthquake_30Days$time[1])
DT <- sapply(Earthquake_30Days$time,FUN=conv.time)
Earthquake_30Days$DateTime <- as.POSIXlt(DT)
#####################################
#Plot Color:Days - Size:Magnitude
#Color Scale
days.from.today <- round(c(Sys.time()-Earthquake_30Days$DateTime)/60,0)
colors.DF <- data.frame(days.from.today,color.scale(days.from.today,color.spec="rgb",extremes=c("red","blue")))
colors.DF <- colors.DF[with(colors.DF, order(colors.DF[,1])), ]
colors.DF$ID <- 1:nrow(colors.DF)
breaks <- seq(1,nrow(colors.DF),length.out=10)
#Size scale
size.DF <- data.frame(Earthquake_30Days$mag,Earthquake_30Days$mag/5)
size.DF <- size.DF[with(size.DF, order(size.DF[,1])), ]
size.DF$ID <- 1:nrow(size.DF)
breaks.size <- seq(0,max(Earthquake_30Days$mag/5),length.out=5)
#Save plot in JPEG
tiff(filename="Earthquake_Map.tif",width=7000,height=4000, res=300)
#Plot
plot(polygons)
plot(Earthquake_30Days, col= colour.scale, cex=Earthquake_30Days$mag/5, pch=16, add=T)
#Title and Legend
title("Earthquakes in the last 30 days",cex.main=3)
legend.pos <- list(x=-28.52392,y=-20.59119)
rect(xleft=legend.pos$x-5, ybottom=legend.pos$y-30, xright=legend.pos$x+30, ytop=legend.pos$y+10, col="white", border=NA)
legendg(legend.pos,legend=c(round(colors.DF[colors.DF$ID %in% round(breaks,0),1],2)),fill=paste(colors.DF[colors.DF$ID %in% round(breaks,0),2]),bty="n",bg=c("white"),y.intersp=0.75,title="Depth",cex=0.8)
text(x=legend.pos$x+5,y=legend.pos$y+5,"Legend:")
legend(x=legend.pos$x+15,y=legend.pos$y,legend=breaks.size[2:5]*5,pch=points(rep(legend.pos$x+15,4),c(legend.pos$y-6,legend.pos$y-9,legend.pos$y-12,legend.pos$y-15),pch=16,cex=breaks.size[2:5]),cex=0.8,bty="n",bg=c("white"),y.intersp=1.1,title="Magnitude")
dev.off()
