GMT

震源地と原子力発電所2

GMT(Generic Mapping Tools) とRのgmt パッケージ

OSはZORIN OS
GMT(Generic Mapping Tools)のバージョンは 4.5.13

インストール用のシェルスクリプト
install_gmt4.sh

netcdf
$ sudo apt-get install libnetcdf-dev

gdal
$ sudo apt-get install gdal-bin
$ sudo apt-get install libgdal-dev

gmt4
$ sudo ./install_gmt4.sh

環境変数の追加: デフォルトではパスが通らない。~/.bashrcの末尾に実行ファイルのパスを追加

(例)インストールしたディレクトリによって異なる。
export GMTHOME=/usr/lib/gmt/
export PATH=/usr/lib/gmt/bin:$PATH

GMT+R+gmt パッケージ 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
library(gmt)
load("in_service.dat") #稼働中の原発  
load("under_construction.dat") #建設中の原発  
load("eq1900_20150121.dat") #地震のデータ   
eqdata <- subset(eq,subset=mag>=6, select=c(latitude,longitude,mag,depth))  
pscoast("-Jm1:25000000 -R120/155/20/50 -W -P -Ba10f5WeSn -Glightgray -K -Y5","quake.eps")
psxy(data.frame(lon=in_service$longitude, lat=in_service$latitude),"-J -R -St0.5 -W0.4 -Gpurple -K -O", "quake.eps")
psxy(data.frame(lon=under_construction$longitude, lat=under_construction$latitude),"-J -R -St0.5 -W0.4 -Gorange -K -O", "quake.eps")
eq0<-subset(eqdata, subset=depth>300)
psxy(data.frame(lon=eq0$longitude, lat=eq0$latitude, mag=10^((eq0$mag-5.9)/2)),"-J -R -Scp -W0.5p,magenta -K -O", "quake.eps")
eq0<-subset(eqdata, subset=depth>50 & depth<=300)
psxy(data.frame(lon=eq0$longitude, lat=eq0$latitude, mag=10^((eq0$mag-5.9)/2)),"-J -R -Scp -W0.5p,blue -K -O", "quake.eps")
eq0<-subset(eqdata, subset=depth>20 & depth<=50)
psxy(data.frame(lon=eq0$longitude, lat=eq0$latitude, mag=10^((eq0$mag-5.9)/2)),"-J -R -Scp -W0.5p,green -K -O", "quake.eps")
eq0<-subset(eqdata, subset=depth<=20)
psxy(data.frame(lon=eq0$longitude, lat=eq0$latitude, mag=10^((eq0$mag-5.9)/2)),"-J -R -Scp -W0.5p,red -K -O", "quake.eps")
gmt.system("pslegend legend.txt -R -J -Dx15.5/0.2/4c/5c/BR -C0.35c/0.35c -F -L2 -G255 -K -O -S","legend3.sh")

legend3.shの実行権限がなしになっていることがある。

ファイルのプロパティ -> パーミッション -> 実行できるようにする

1
2
3
4
5
system("./legend3.sh>>quake.eps")
system("psscale -D8/0/15c/0.5ch -Cdepth.cpt -Ba100f50g100 -K -O -Y-1>>quake.eps")
pstext(data.frame(Lon=150,Lat=33.5,Size=15,Angle=0,Font=21,Justify="CM",Text="Magnitude"),"-J -R -K -O", "quake.eps")
pstext(data.frame(Lon=137.5,Lat=16,Size=15,Angle=0,Font=21,Justify="CM",Text="Depth"),"-J -R -N -K -O", "quake.eps")
pstext(data.frame(Lon=137.5,Lat=53,Size=17,Angle=0,Font=21,Justify="CM",Text="Earthquake M>=6.0 (1900_2015/1/21) & Nuclear power stations"),"-J -R -N -O", "quake.eps")

gmt.system関数が使えなくてもsystemがある。

legend.txtの内容

N 1
S 0.5c c 1.12p 255/0/0 0.5p 1.5c mag6.0
S 0.5c c 3.55p 255/0/0 0.5p 1.5c mag7.0
S 0.5c c 11.22p 255/0/0 0.5p 1.5c mag8.0
S 0.5c c 35.48p 255/0/0 0.5p 1.5c mag9.0

depth.cptの内容

0 255 0 0 20 255 0 0
20 0 255 0 50 0 255 0
50 0 0 255 300 0 0 255
300 255 0 255 720 255 0 255
B 0 0 0
F 255 255 255
N 128 128 128

Inkscapeでpng ファイルに変換