コンピュータを研究に使うための私的メモ。Python、Fortran、Ubuntu、etc...


基本

インポート

import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

m = Basemap(...)                   # 地図を作成(指定しない場合はax=plt.gca()に描かれる)
m.drawmapboundary()           # 外枠線を描く
m.drawcoastlines()                 # 海岸線を描く

m.drawmeridians(...)               # 経度線とラベル
m.drawparallels(...)                 # 緯度線とラベル

# 描画関係は基本的にmatplotlibと共通
m.contour(...)
m.contourf(...)
m.plot(...)

プロット

作成したBasemapクラスのインスタンスに緯度経度の配列を渡すことで図の中の座標に変換される。
m = Basemap(projection='merc', llcrnrlat=-90, urcrnrlat=90, \
                llcrnrlon=0, urcrnrlon=360, resolution='c')
X, Y = meshgrid(lon,lat)
x, y = m(X, Y)                    #projection='cycl' の場合はこの操作はなくてもよい。

m.contourf(x, y, z)
projection='cycl' で普段プロットしていると投影法を変えた時に何もプロットされないことがある。 m(X, Y) による変換を忘れていることが多いので気を付けること。

地図を描く

map = Basemap(keyward**)
キーワードで地図の種類、範囲等を指定して地図を作成する。

地図の種類

projection地図の投影法
以下の他にも面積や距離が等しくなるとかいろいろある。こちらを参考。
  • 正距円筒図法(cyl
  • メルカトル図法("merc")
|
  • 投影図法(stere)
  • 極投影図法(npstere)
  • GrADSとの対応
set gxout latloncyl
set gxout npsnpstere/spstere

範囲の指定

  • 両端の値で指定する場合
llcrnrlon左下端の経度(degree)
urcrnrlon右上橋の経度
llcrnrlat左下端の緯度
urcrnrlat右上端の緯度
  • 幅で指定する場合
width地図の東西幅(メートル)
heiht地図の南北幅
lon_0地図の中心経度(degree)
lat_0地図の中心緯度
  • 極中心の場合(npstere,spstere,nplaea,splaea,npaeqd,spaeqd)
lon_0中心となる経度(GrADSでは180度)
boundinglat境界の緯度
地図の種類によっては内部で計算されるので指定しなくてよいものもある。

その他

resolution地図の境界の解像度。c(crude) < l(low) < i(intermediate) < h(high) < f(full) の順で細かくなる。デフォルトはc
aspect地図のアスペクト比。(y方向, x方向)
fix_aspect地図の縦横比を地形に対して固定するか。デフォルトはTrue
lat_tsスケールが正しくなる緯度。stereographicかmercatorの場合に指定できる。stereographicではlat_0、mercatorでは0がデフォルト
round極中心の図のときに枠を円にするか。

境界の設定

Basemapのインスタンスにメソッドで指定していく。

境界の描画

#海岸線を描く
map.drawcoastlines()
#国境を描く
map.drawcountries()
  • キーワード
linewidth線幅
color
antialiasedアンチエイリアスをかけるか。デフォルトはTrue

大陸の色

map.fillcontinents(color='coral', lake_color='aqua')
  • キーワード
color大陸の色
lake_color湖の色

背景色

背景(海)の塗りつぶしと地図の枠を引く
map.drawmapboundary(fill_color='aqua') 
  • キーワード
linewidth枠線の幅
color枠線の色
fill_color塗りつぶす色

緯度・経度線の設定

渡す緯度、経度の配列はデータの配列やBasemapで指定する範囲と一致しなくてもいい。0-360度、-90-90度の配列を使ってあとは間隔だけを考えればいい。
#経度線を引く(0度から30度間隔)
map.drawmeridians(np.arange(0, 360, 30), labels=[True, Flase, False, True])
#緯度線を引く(-90度から30度間隔)
map.drawparallels(np.arange(-90, 90.1, 10), labels=[False, True, True, False])
  • キーワード
color線色
linewidth線幅
dashed点線のパターン。デフォルトは[1, 1]。つまり1ピクセル引いて1ピクセル引かないを交互に繰り返す。
labels4要素のリストで指定。それぞれ左枠、右枠、上枠、下枠との交点にラベルをつけるかを指定

バグ

(2012.5.21) basemap-1.0.2(以下?)では緯度線の指定に整数以外の配列を渡すとエラーになるバグがあるらしい。1.0.3では修正されている。
http://www.mail-archive.com/matplotlib-users@lists...

緯度・経度線を引かずに緯度・経度のラベルだけ設定する。

drawmeridiansやdrawparallelsは、メルカトル図法などの等間隔でない地図の場合でも特定の緯度・経度にラベルを引けるのでとても便利だが、緯度・経度のグリッドが図に入ってしまうのがときどき邪魔に感じる時がある。GrADSや通常のmatplotlibのようにグリッドではなく軸メモリを入れるには次のようにFormatterを作って設定する。

こちらを参考。緯度・経度の値(130とか)を文字列(130Eとか)に変換する関数(lon2txtはBasemapのコードを参考にして作った。

import matplotlib.ticker
class BasemapYaxisFormatter(matplotlib.ticker.Formatter):
    def __init__(self,baseMap):
        self.baseMap=baseMap
    def __call__(self,y,pos=1):
        lon,lat = self.baseMap(0,y,inverse=True)
        return lat2txt(lat,pos=pos)

class BasemapXaxisFormatter(matplotlib.ticker.Formatter):
    def __init__(self,baseMap):
        self.baseMap=baseMap
    def __call__(self,x,pos=1):
        lon,lat = self.baseMap(x,0,inverse=True)
        return lon2txt(lon,pos=pos)
このFormatterをaxesに設定する。
m=Basemap(......., suppress_ticks=False)    #suppress_ticksをFalseにして軸メモリをオンにする。
ax=gca() 
ax.xaxis.set_major_formatter(BasemapXaxisFormatter(m))
ax.yaxis.set_major_formatter(BasemapYaxisFormatter(m))

その他

大円を引く

map.drawgreatcircle(lon1, lat1, lon2, lat2)

プロット

matplotlibの2次元プロットとほぼ同じ。まず位置の座標を表すnumpy2次元配列をつくり、それをBasemapインスタンスの引数に与えて地図における2次元グリッド配列のようなものをつくる。
m = Basemap(projection='cyl', llcrnrlat=-90, urcrnrlat=90, \
                llcrnrlon=0, urcrnrlon=360, resolution='c')

lon = np.arange(-180.,181.,60.)
lat = np.arange(-90.,91.,30.)

#2次元グリッド配列をつくる。計算はコレを利用する。
lons, lats = np.meshgrid(lon, lat)
Z = function(lons, lats)

#描画したい2次元配列の座標に対応する2次元配列の地図での座標を求める。
#座標配列は緯度、経度で表しておく。
X, Y = map(lons, lats)

map.countour(X, Y, Z)
プロットの線種、間隔等のキーワードはmatplotlibと共通。

コンター

map.contour(X, Y, Z)
コンターのラベル
コンターのラベルはplt.clabelで設定できるが配列の範囲とmapの範囲が合っていないと図の外側にまでラベルがはみ出してしまう。これは図の範囲の情報をcontourが知らないかららしい。mask配列を使って解決する(参考。だいぶましになるけど完全には解決しない。アップデートに期待。バージョン1.0.6で図の外側はマスクするように改善されたようだ。
m = Basemap(...)

X, Y = m(lon, lat)
mask = X < m.xmin | X > m.xmax | Y > m.ymax | Y < m.ymin
z = np.ma.array(z,mask=mask)
CS= m.contour(X,Y,z,np.hstack([np.arange(-2,0,0.02),np.arange(0.02,2,0.02)]),colors='k')

plt.clabel(CS,fmt='%g')
コンターを一定間隔で引く
よく使うのでメモ
from pylab import *
#100間隔でコンターを引く
m.comtour(X, Y, var, locator=MultipleLocator(100))  

シェード

from pylab import *
map.contourf(X, Y, Z)

#カラーマップを指定して100間隔で描画
map.contourf(X, Y, Z, cmap=cm.RdBu, locator=MultipleLocator(100))

カラー

map.imshow(Z)
orginの位置とかは自動的に指定されるらしい。

軌跡

xtrack = [100, 200, 300, .....]
ytrack = [50, 20, 100, ....]
map.plot(xtrack, ytrack, '-o')

緯度経度をサイクリックな値で指定する。

緯度、経度の情報はプロットする関数(plot, contour, contourfなど)に1番目、2番目の変数として与える。この値がマイナスで入っているとBasemapで指定したマップの範囲がプロットの際にサイクリックに解釈されない。
例えば以下のように130Eから110Wの範囲を指定したBasemapオブジェクトに西経をマイナスで表したlonを使ってプロットしてみる。
from pylab import *

#変数varには-180度から180度まで東から西にデータが入っている。
lon = np.arange(-180,180,1.25)
lat = np.arange(-90,90,1.25)

m = Basemap(projection='cyl', llcrnrlat=0, urcrnrlat=60, \
           llcrnrlon=130, urcrnrlon=250, resolution='c',lon_0=180)
m.drawcoastlines(linewidth=0.5)
m.drawmapboundary()

m.drawmeridians(np.arange(0,360,30),labels=[0,0,0,1],linewidth=0.5,fontsize=14)
m.drawparallels(np.arange(-90,90,10),labels=[1,0,0,0],linewidth=0.5,fontsize=14)

x, y = np.meshgrid(lon, lat)
X, Y = m(x, y)

CF=m.contourf(X, Y, var[:,:]) 
CB=m.colorbar(CF)
すると180以上の値がlonに入っていないためなのか西半球がプロットされない(マップでは正しく解釈されて地図は描かれている)。



Basemapのところで西半球をマイナスの緯度で表現して、
m = Basemap(projection='cyl', llcrnrlat=0, urcrnrlat=60, \
           llcrnrlon=130, urcrnrlon=-110, resolution='c',lon_0=180)
とするとマイナスで指定した側が図の左側にきて110W-130Eで描かれてしまう。



対策(ver. 1.0.5以上)
Basemap 1.0.5からは緯度経度がプロットの際にサイクリックに解釈されるようにプロットする関数にオプションlonlatが追加されている。
m = Basemap(projection='cyl', llcrnrlat=0, urcrnrlat=60, \
           llcrnrlon=130, urcrnrlon=250, resolution='c',lon_0=180)
m.drawcoastlines(linewidth=0.5)
m.drawmapboundary()

m.drawmeridians(np.arange(0,360,30),labels=[0,0,0,1],linewidth=0.5,fontsize=14)
m.drawparallels(np.arange(-90,90,10),labels=[1,0,0,0],linewidth=0.5,fontsize=14)

x, y = np.meshgrid(lon, lat)
X, Y = m(x, y)

CF=m.contourf(X, Y, var[:,:], latlon=True) 
とすると望み通りサイクリックに解釈されて図がかける。

カラーバー

カラーバーを描く

プロットの返り値を引数にしてカラーバーを描画する。
CF = map.contourf(X, Y, Z)
map.colorbar(CF, location='right', size='5%')
  • キーワード
locationカラーバーの位置。top、bottom、left、right
sizeカラーバーの幅。地図の幅(高さ)に対して整数のパーセントで指定。
padカラーバーと地図との距離。地図の幅(高さ)に対して整数パーセントで指定。
ticksカラーバーの目盛り。リストで。

カラーバーのラベル

colorbarメソッドはmatplotlibのcolorbarインスタンスを返すのでmatplotlibと同じようにする。
CB = map.colorbar(CF)
CB.set_clabel('label')

海や陸をマスクする

maskoceans()という関数を使うと海部分をmaskしたmasked arrayを求めることができる。
from mpl_toolkits.basemap import maskoceans
#西経はマイナスで表さないと認識してくれないので変換
lon[lon>180] = lon[lon>180] - 360
var = maskoceans(lon, lat, datain)
海をマスクするにはmaskを入れ替える。
from mpl_toolkits.basemap import maskoceans

lon[lon>180] = lon[lon>180] - 360
var = maskoceans(lon, lat, datain)
var.mask = ~var.mask

グリッドに内挿する

interpという関数で緯度経度のグリッドデータを別の緯度経度のグリッドに内挿することができる。
from mpl_toolkits.basemap import interp

#datain                  :  内挿する2次元のarray datain.shape = (len(lat), len(lon))
#lonin, latin         :  内挿前のデータ(datain)の緯度、経度(1次元array)
#lonout, latout  :  内挿するグリッド。それぞれ2次元配列(meshgrid等でつくる)
data = interp(datain, lonin, latin, lonout, latout)
  • オプション
order内挿方法。nearest-neighbor(order=0)、bilinear(order=1)、cubic spline(order=3)。デフォルトはorder=1
×

この広告は60日間更新がないwikiに表示されております。

情報

管理人のみ編集できます