Article Outline

原発と地震 (その2) 日本と台湾

使用したデータ

マグニチュード５以上の地震データはアメリカ地質調査所(USGS)

Rコード

パッケージの読み込み。関数読み込み。など

``````library(oce)
data(coastlineWorldFine, package="ocedata")
#
#
# stackoverflow : Drawing a Circle with a Radius of a Defined Distance in a Map
#
mapLonLat <- function(LonDec, LatDec, Km) {
ER <- 6371 #Mean Earth radius in kilometers. Change this to 3959 and you will have your function working in miles.
AngDeg <- seq(1:360) #angles in degrees
Lat2Deg <- Lat2Rad*(180/pi)#Latitude of each point of the circle rearding to angle in degrees (conversion of radians to degrees deg = rad*(180/pi) )
Lon2Deg <- Lon2Rad*(180/pi)#Longitude of each point of the circle rearding to angle in degrees (conversion of radians to degrees deg = rad*(180/pi) )
return(data.frame(Lon2Deg,Lat2Deg))
}
#
nuclearPower<- nuclearPower[order(nuclearPower\$Country),]
# In service と Under construction　の原発のみ(Decommissioned以外)
nuclearPower<- nuclearPower[nuclearPower\$CurrentStatus != "Decommissioned"  ,]``````

原発別（１６６）pngファイル作成

``````for (i in 1:nrow(nuclearPower)){
np<- c(nuclearPower\$longitude[i], nuclearPower\$latitude[i])
LonDec<- np[1]
LatDec<- np[2]
Km<- 120
LonLat<- mapLonLat(LonDec, LatDec, Km)
par(mar=c(2, 2, 3, 2))
lonlim <- range(LonLat[,1])
latlim <- range(LonLat[,2])
#
# 正距方位図法  azimuthal equidistant projection
aeqd_proj <- paste("+proj=aeqd +lon_0=",np[1]," +lat_0=",np[2])
#
png(paste0("npeq",i,".png"),width=800,height=800)
mapPlot(coastlineWorldFine, projection=aeqd_proj ,
col="lightgray", longitudelim=lonlim, latitudelim=latlim)
title(paste(nuclearPower[,1][i],"(",nuclearPower[,2][i],")"))
mapPoints(np[1],np[2],pch=24,bg="darkgreen",col="black",cex=1.5)
#
angle<- 1
pos<- 3
Km<- 30
dist<- mapLonLat(np[1],np[2],Km)
mapPolygon(dist[,1],dist[,2],border="blue", lwd=0.5)
mapText(dist[angle,1],dist[angle,2],labels=paste(Km ,"km"), cex =1,col="blue",pos=pos)
Km<- 60
dist<- mapLonLat(np[1],np[2],Km)
mapPolygon(dist[,1],dist[,2],border="blue", lwd=0.5)
mapText(dist[angle,1],dist[angle,2],labels=paste(Km ,"km"), cex =1,col="blue",pos=pos)
Km<-100
dist<- mapLonLat(np[1],np[2],Km)
mapPolygon(dist[,1],dist[,2],border="blue", lwd=0.5)
mapText(dist[angle,1],dist[angle,2],labels=paste(Km ,"km"), cex =1,col="blue",pos=pos)
#
mapPoints(eqdata\$longitude,eqdata\$latitude,pch=21,bg=rgb(1,0,0,alpha=1),col="gray20",cex=1.2)
dev.off()
}``````