Sys.time()

[1] “2014-10-27 12:48:44 JST”

data.location<-2 #1:local 2:web
data.file<-"all_month.csv"
if(data.location==1){
username<-Sys.info()['user']
path01<-paste("C:/Users/",username,"/Desktop/Earthquake_Data/",sep="")
dataset<-read.table(file=paste(path01,data.file,sep=""),sep=",",header=T,as.is=T,skip=0)
}else{
dataset<-read.csv(
paste("http://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/",data.file,sep=""),
header=T,
as.is=T,
skip=0
) 
}
dataset.j<-dataset[grep("Japan",dataset[,14]),]
#Prepare for googleVis
dataset.j$latlong<-paste(dataset.j$latitude,dataset.j$longitude,sep=":")
dataset.j$jst<-as.POSIXct(sub("Z","",sub("T"," ",dataset.j$time)))+3600*9 #UTC+9
dataset.j$info<-paste(
"Time(UTC)=",dataset.j$time,
", UTC+9=",dataset.j$jst,
", Depth(km)=",dataset.j$depth,
", Mag=",dataset.j$mag,
", MagType=",dataset.j$magType,
", Place=",dataset.j$place
,sep="")
dataset.j$energy.GJ<-round(10^((dataset.j$mag)*1.5+4.8)*(10^-9),4)
gMap<-gvisMap(
dataset.j,"latlong","info",
options=list(
showTip=TRUE, 
showLine=TRUE, 
enableScrollWheel=TRUE,
mapType='terrain', #'hybrid','normal','terrain','satellite'
useMapTypeControl=TRUE
)
)
#Prepare for maps
for(mmm in 1:length(dataset.j[,1])){
tmp<-dataset.j$mag[mmm]+10^-3
if(round(tmp)==floor(tmp)){tmp<-floor(tmp)}else{tmp<-(round(tmp)-0.5)} 
dataset.j$magMap[mmm]<-tmp
}
japan<-map(xlim=c(120,150),ylim=c(20,50))
# http://cran.r-project.org/web/packages/maps/maps.pdf#page=8
#str(japan)
japan<-as.data.frame(japan[c("x","y")])
japan.map<-ggplot(japan,aes(x,y))
japan.map<-japan.map+geom_path()
# http://cran.r-project.org/web/packages/ggplot2/ggplot2.pdf#page=79
japan.map<-japan.map+geom_point(
data=dataset.j,
aes(
x=longitude,
y=latitude,
size=magMap
)
)
japan.map<-japan.map+scale_size_continuous(paste("Magnitude(M)","\n","M~(M+0.5)"))

Earthquake Map,maps

plot(japan.map)

Earthquake Map,google map

print(gMap,tag="chart")
#prepare for 3D-plot
lon.range<-c(120,150)
lat.range<-c(20,60)
Latitude.degree<-dataset.j$latitude
Longitude.degree<-dataset.j$longitude
Magnitude<-dataset.j$mag
ang<-c(110) # angle
z.border<-4.5
z.min<-max(floor(min(Magnitude)),z.border)
z.max<-max(Magnitude)
ddd<-1
map2d<-map(xlim=c(lon.range[1],lon.range[2]),ylim=c(lat.range[1],lat.range[2]))
map2d<-as.data.frame(map2d[c("x","y")])
map2d<-cbind(map2d,z=z.min)

Earthquake Map,3D-Plot

ppp<-1
for(ccc in 1:length(map2d[,1])){
if(is.na(map2d[,1][ccc])==T){
point.y=map2d$y[ppp:ccc-1]
point.x=map2d$x[ppp:ccc-1]
point.z=map2d$z[ppp:ccc-1]
map.line<-scatterplot3d(
y=point.y,
x=point.x,
z=point.z,
angle=ang[ddd],
lwd="2",
type="l",
scale.y=1,
color=2,
xlim=c(lon.range[1],lon.range[2]),
ylim=c(lat.range[1],lat.range[2]),
zlim=c(z.min,z.max),
xlab="",
ylab="",
zlab="",
x.ticklabs="",
y.ticklabs="",
z.ticklabs=""
)
ppp<-ccc
par(new=T)
gc()
gc()
}
}
earthquake.point<-scatterplot3d(
y=Latitude.degree,
x=Longitude.degree,
z=Magnitude,
angle=ang[ddd],
type="h",
lwd="3",
scale.y=1,
color=4,
main=paste("In the past 30days-Magnitude-over",z.border,sep=""),
xlim=c(lon.range[1],lon.range[2]),
ylim=c(lat.range[1],lat.range[2]),
zlim=c(z.min,z.max),pch=" "
)

Histogram

PlotTitle<-"Japan"
bw=c(5,5,10,0.1)
for(ppp in 1:4){ 
switch(ppp,
g<-ggplot(dataset.j,aes(x=latitude)),
g<-ggplot(dataset.j,aes(x=longitude)),
g<-ggplot(dataset.j,aes(x=depth)),
g<-ggplot(dataset.j,aes(x=mag))
)  
g<-g+geom_histogram(alpha=0.2,binwidth=bw[ppp],colour="blue")
g<-g+labs(title=PlotTitle)
g<-g+theme(plot.title=element_text(size=9))
g<-g+theme(axis.title.x=element_text(size=9),axis.title.y=element_text(size=9)) 
g<-g+theme(axis.text.x=element_text(size=9),axis.text.y=element_text(size=9)) 
switch(ppp,
lat<-g,
lon<-g,
dep<-g,
magn<-g
)
}
grid.arrange(lat,lon,dep,magn,nrow=1,ncol=4)

## [1] "Time(UTC)=2014-10-11T02:35:46.220Z, UTC+9=2014-10-11 11:35:46, Depth(km)=13.48, Mag=6.3, MagType=mwp, Place=154km ENE of Hachinohe, Japan"

Earthquake List in the past 30 days

E.list<-dataset.j[,c(17,5,19,6,14,2,3,4)]
E.list<-data.frame(No=c(1:nrow(E.list)),E.list)
gTable<-gvisTable(E.list)
print(gTable,tag="chart")

Total Energy

totalEnergy<-sum(dataset.j[,19])
cat("Total Energy(TJ)=",totalEnergy*10^-3,"\n")
## Total Energy(TJ)= 271.1198
cat("Energy of ",round(totalEnergy*(10^-3)/67,1),"Little Boy")
## Energy of  4 Little Boy

Covariate Plot

Data Source/Reference