R

2015日本周辺震源マップ

marmap ,RCurl パッケージ

marmap パッケージについては
海底地形と地形データでも記事にしています。

震源のデータはUSGSよりRCurl パッケージを用いて入手します。

（参考）
API Documentation - Earthquake Catalog

2015日本周辺震源マップ

コードは３つの部分に分かれています。

  1. 地形データをNOAAからダウンロード(getNOAA.bathy 関数)
  2. 震源データダウンロード(getURL 関数)。保存。
  3. 震源マップ作成

一括処理（震源マップ作成）

地形データをNOAAからダウンロード

1
2
3
4
5
6
7
8
9
10
11
library(marmap)
library(RCurl)
library(colorspace)
lon<-c(119,155)
lat<-c(19,51)
papoue <- getNOAA.bathy(lon1 = lon[1], lon2 = lon[2],lat1 = lat[1], lat2 = lat[2], resolution = 8)
oce<-c(rgb(0,0,0),rgb(0,5/255,25/255),rgb(0,10/255,50/255),rgb(0,80/255,125/255),
rgb(0,150/255,200/255),rgb(86/255,197/255,184/255),rgb(172/255,245/255,168/255),
rgb(211/255,250/255,211/255),rgb(250/255,1,1)) 
#land<-terrain.colors(12)
land<-terrain_hcl(12)

一括処理(USGS:地震データダウンロード。保存。)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
#震源データをダウンロードする年（ここでは2015年のみ）
start<-2015 ; end<-2015
library(RCurl)
for (i in start:end){
year<-as.character(i)
#format=csvとする
url<-paste0("http://earthquake.usgs.gov/fdsnws/event/1/query?format=csv&starttime=",as.numeric(year),"-01-01&endtime=",as.numeric(year)+1,
"-01-01&minmagnitude=4&minlatitude=",lat[1],"&maxlatitude=",lat[2],"&minlongitude=",lon[1],"&maxlongitude=",lon[2],"&eventtype=earthquake")
doc <- getURL(url)
filename<-paste0("USGS",year,"_",lon[1],"_",lon[2],"_",lat[1],"_",lat[2],".csv")
out <- file(filename, "w") # 書き込みモードで開く
writeLines(doc,out,sep = "\n") 
close(out)   
}

一括処理（震源マップ作成）

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
#start<-2015 ; end<-2015
#lon<-c(119,155)
#lat<-c(19,51)
for (i in start:end){
year<-as.character(i)
filename<-paste0("USGS",year,"_",lon[1],"_",lon[2],"_",lat[1],"_",lat[2],".csv")
eq<-read.csv(filename)
#eq<-eq[paste(eq$type)=="earthquake",]
eqdata <- subset(eq, select=c(time,latitude,longitude,depth,mag,nst,gap,dmin,rms))
#summary(eqdata)
#
#eqdata<- subset(eq, grepl(year,eq$time))
#マグニチュードの昇順に並べ替え
sortlist <- order(eqdata$mag)
eqdata<- eqdata[sortlist,]
#Depthによって色を変えるための準備
for (i in 1:nrow(eqdata)){
if (eqdata$depth[i]<=20)  { 
eqdata$col[i] <- rgb(1,0,0,alpha=0.4) 
} else if ((eqdata$depth[i]>20) && (eqdata$depth[i]<=50)) {
eqdata$col[i] <- rgb(0,1,0,alpha=0.4)
} else if ((eqdata$depth[i]>50) && (eqdata$depth[i]<=300)) {
eqdata$col[i] <- rgb(0,0,1, alpha=0.4)
} else if (eqdata$depth[i]>300) {
eqdata$col[i] <- rgb(1,0,1, alpha=0.4)
}
}
#### 震源マップ保存
png(paste0("USGS",year,"_",lon[1],"_",lon[2],"_",lat[1],"_",lat[2],".png"),width=1000,height=800)
#
par(mar=c(5, 4, 4, 3), xpd=TRUE)
plot(papoue, image = TRUE, land = TRUE,axes=F,xlab="",ylab="", lwd = 0.1,
bpal = list(c(0, max(papoue), land),c(-8000,0,oce)))
# Making the coastline more visible
plot(papoue, deep = 0, shallow = 0, step = 0,lwd = 0.4, add = TRUE)
#震源地をプロット
magsize<-(eqdata$mag-3)
points(eqdata$longitude,eqdata$latitude,pch=21,cex =magsize, col="black",bg=eqdata$col)  
cexsize<-(c(4,5,6,7,8,9)-3)
#凡例（サイズ：マグニチュード）
legend(x=lon[2]+1,y=(lat[1]+lat[2])/2+12,title="Magnitude",legend = c("4","5","6","7","8","9"),  
cex=1,pch=21,col="black",pt.bg=rgb(1,0,0,alpha=0.6),pt.cex =cexsize,bg ="gray", x.intersp = 2, y.intersp = 2.5)
#凡例（色：震源の深さ）
legend(x=lon[2]+1,y=(lat[1]+lat[2])/2-5,title="Depth",legend = c("  0～20 km","20～50 km","50～300km","300km～"),  
cex=1,pch=21,col="black",pt.bg=c(rgb(1,0,0,alpha=0.6),rgb(0,1,0,alpha=0.6),rgb(0,0,1,alpha=0.6),rgb(1,0,1,alpha=0.6)),
pt.cex =3.2,bg ="gray", x.intersp = 1.5, y.intersp = 2)
title(paste("Earthquakes (",year,")"),cex.main=2)
#軸を描く
segments(lon[1], lat[1], lon[1], lat[2],lwd=2)
segments(lon[1], lat[1], lon[2], lat[1],lwd=2)
for (i in 1:(floor((lat[2]-lat[1])/10)+1)){
	segments(lon[1], 10*(ceiling(lat[1]/10)+(i-1)), lon[1]-0.5, 10*(ceiling(lat[1]/10)+(i-1)),lwd=0.5)
	text(lon[1]-1.5, 10*(ceiling(lat[1]/10)+(i-1)),as.character(10*(ceiling(lat[1]/10)+(i-1))))
 }
for (i in 1:(floor((lon[2]-lon[1])/10)+1)){
	segments(10*(ceiling(lon[1]/10)+(i-1)),lat[1]-0.5, 10*(ceiling(lon[1]/10)+(i-1)),lat[1],lwd=0.5)
	text(10*(ceiling(lon[1]/10)+(i-1)),lat[1]-1.2, as.character(10*(ceiling(lon[1]/10)+(i-1))))
 }
text((lon[1]+lon[2])/2,lat[1]-2.2,"Longitude")
text(lon[1]-2.7,(lat[1]+lat[2])/2,"Latitude",srt=90)
dev.off()
}