memoring

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 latlon	cyl
set gxout nps	npstere/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ピクセル引かないを交互に繰り返す。
labels	4要素のリストで指定。それぞれ左枠、右枠、上枠、下枠との交点にラベルをつけるかを指定

(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の２次元プロットとほぼ同じ。まず位置の座標を表すnumpy２次元配列をつくり、それを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次元配列の座標に対応する２次元配列の地図での座標を求める。
#座標配列は緯度、経度で表しておく。
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