(Tips) 地図上に点をプロット

作成者:山崎哲

概要

地図上に点をプロットするスクリプト。メルカトル図法用とポーラーステレオ用の2つを紹介する。

使い方

プロットしたい経度・緯度情報の入ったテキストファイルを用意する。中身は、

%less points.txt
 116.26 39.55
 77.11 28.33
 139.46 35.40
 126.58 37.32

というフォーマットで、経度、緯度の順に1スペース開けて並べる。このフォーマットは、fortranで

do i = 1, imax
write(1,*) lon(i), lat(i)
enddo

と出力したformattedファイルと同じ。imaxはプロットする点の数を示す。

テキストファイルへのパスを記述して、rubyスクリプト(以下サンプル)を実行する。

サンプル

メルカトル図法版

require "numru/ggraph"
require "numru/dcl"
include NumRu

## setting

# file
file1 = File.open("points.txt",'r') #読み込むテキストファイル名

# IO switch
io = 2 #1 standard, 2 ps, 4+[d] png

# draw range
x1 = 0
x2 = 360
y1 = -90  # y2>y1
y2 = 90

# label interval
lx = (x2-x1)/4
ly = (x2-x1)/8

# mark
marksize=0.01
marktype=2
markindex=20

# general
DCL.gropn(io)
DCL::grfrm

# mapping
DCL::grswnd(x1, x2, y1, y2) #描画範囲を指定
DCL::grsvpt(0.2, 0.9, 0.2, 0.7) #フレームの大きさを指定
DCL::grstrn(1)
DCL::grstrf

#x-label
DCL::uxaxdv('B', lx, lx)
DCL::uxaxdv('T', lx, lx)
DCL::uxsttl('B', 'Lon.', 0.0)

# y-label
DCL::uyaxdv('L', ly, ly)
DCL::uyaxdv('R', ly, ly)
DCL::uysttl('L', 'Lat.', 0.0)
DCL::uddclv(0.0)

# draw map & title
DCL::sglset('lclip','true') 
DCL::umpmap('coast_world')
DCL::uxmttl('T', 'Map & points', 0.0)

# draw mark
DCL::uusmks(marksize)
DCL::uusmkt(marktype)
DCL::uusmki(markindex)
loop do
p0 = file1.gets
if p0.nil? then # stop if file reaches EOF
break
end

up = p0.split(/\s+/)
upx = up[1].to_f
upy = up[2].to_f
 print upx, ', ', upy, ' '
xmark = NArray.sfloat(1)
ymark = NArray.sfloat(1)
xmark[0] = upx
if upx < 0 then
xmark[0] = upx + 360
end
ymark[0] = upy
DCL.uumrk(xmark,ymark) 
end

DCL.grcls

plot_Mercator.png

ポーラーステレオ版

require "numru/ggraph"
require "numru/dcl"
include NumRu
include GGraph

## setting
fname='points.txt' # テキストファイルの名前
fname_title = fname.split(/_/)[0]
title='Map & points' #title 
# IO switch
io = 1 #1 standard, 2 ps, 4+[d] png
# draw range
x1 = -180
x2 = 180
y1 = 30 #y2>y1
y2 = 90

# viewport
view_r = (y1-y2).abs  
view_center_lon = (x1+x2)/2
view_center_lat = y2
if y1.abs > y2.abs then
view_center_lat = y1
end
# label
label_B = view_center_lon.to_s
if(y1 < 0) then
label_B = (view_center_lon+180).to_s
end
#latlon line inverval
dmj = (x2-x1)/8.0  # Major lineの指定。経度方向に8分割。
dmn = dmj/6.0 # Minor lineの指定。Major lineを6分割。
# mark
marksize=0.01
marktype=2
markindex=20

# file
file1 = File.open(fname,'r')

# general
DCL.gropn(io)
DCL::grfrm
# mapping
DCL::grsvpt(0.2, 0.9, 0.2, 0.7)
DCL::grstrn(31)
DCL::grswnd(x1, x2, y1, y2)
DCL::umscwd(view_center_lon, view_center_lat, view_r) # (lon, lat, r)円形のウィンドウを設定
DCL::umpfit
DCL::grstrf
DCL::slpvpr(1)
# draw map & title
DCL::sglset('lclip','true') 
DCL::umpmap('coast_world')
DCL::umpset('dgridmj',dmj) 
DCL::umpset('dgridmn',dmn) 
DCL::umpset('dgrplmn',10.0)       # Minor lineを極から何度離れた緯度から描くか。
DCL::umpglb                       # 緯経線と地球の輪郭を描く

# draw mark
DCL::uusmks(marksize)
DCL::uusmkt(marktype)
DCL::uusmki(markindex)
loop do
p0 = file1.gets
if p0.nil? then # stop if file reaches EOF
break
end
up = p0.split(/\s+/)
upx = up[1].to_f
upy = up[2].to_f
 print upx, ', ', upy, ' '
xmark = NArray.sfloat(1)
ymark = NArray.sfloat(1)
xmark[0] = upx
if upx < 0 then
xmark[0] = upx + 360
end
ymark[0] = upy
if upy>=y1 && upy<=y2 then 
DCL.uumrk(xmark,ymark) 
end
end
#title
if y1 > 0 then
 11.times do
 DCL::uxsttl('B', ' ', 0)  #タイトルを描く位置を調整。
 end
else
 5.times do
 DCL::uxsttl('B', ' ', 0)
 end
end
DCL::uxsttl('B', label_B+'E', 0)
if y1 > 0 then
 4.times do
 DCL::uxmttl('T', ' ', 0)
 end
else
 8.times do
 DCL::uxmttl('T', ' ', 0)
 end
end
DCL::uxsttl('T', fname_title, -1)
DCL::uxmttl('T', title, 0)

DCL.grcls

plot_NP.png

南半球については,

y1 = -90
y2 = -40

などとする:

plot_SP.png

メモ

  • テキストファイルは実数でも整数でもOK。ゼロが先頭にくる(-001.23)場合もOK。
  • 経度は東経だけで書いても(220.)、西経をマイナスで記述しても(-140.)OK。
更新日時:2013/04/29 11:22:28
キーワード:[プロット] [サンプル]
参照: