15
21

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

Pythonで世界地図(matplotlib basemap toolkit翻訳)-4

Last updated at Posted at 2018-07-24

Plotting data on a map (Example Gallery)

注意)
データのURLにアクセスしてもなかったり、データが新しくなりファイル名が変わっている可能性があります。

次に、Basemapインスタンスメソッドを使用してデータをマップにプロットする方法を示す一連の例を示します。より多くの例は、ベースマップソース配布のexamplesディレクトリに含まれています。データをプロットするためのいくつかのBasemapインスタンスメソッドがあります:

・ contour(): 等高線を描きます。
・ contourf(): 塗りつぶされた輪郭を描く
・ imshow(): 画像を描画する。
・ pcolor(): 疑似カラープロットを描きます。
・ pcolormesh(): pseudocolorプロット(通常のメッシュの場合はより高速版)を描画します。
・ plot(): 線および/またはマーカーを描画する。
・ scatter(): マーカーで点を描く。
・ quiver(): 描画ベクトル。
・ barbs(): 風の棘を引く。(風向・風速の記号)
・ drawgreatcircle(): 大きな円を描きなさい。(大圏コース)

これらのインスタンスメソッドの多くは、対応するmatplotlib Axesインスタンスメソッドに単純に転送し、いくつかの追加のプリ/ポスト処理と引数のチェックを行います。また、ベースマップに関連付けられたAxesインスタンスを使用して、matplotlib pyplotインターフェイスまたはOO APIを使用してマップに直接プロットすることもできます。

Basemapインスタンスメソッドの使い方の詳細については、「Matplotlib Basemap Toolkit API」を参照してください。

以下に例を示します(その多くは、httpを介してデータセットを取得するためにnetcdf4-pythonモジュールを利用しています)。

ベースマップ上に等高線を描画する
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
# 50N、100Wを見下ろす衛星の視点を用いて正射投影図投影を設定します。
# 低解像度の海岸線を使用する。

map = Basemap(projection='ortho',lat_0=45,lon_0=-100,resolution='l')
# 海岸線を描く、国境を塗る、大陸を塗る
map.drawcoastlines(linewidth=0.25)
map.drawcountries(linewidth=0.25)
map.fillcontinents(color='coral',lake_color='aqua')
# 地図投影領域(投影肢)の端を描画し、
map.drawmapboundary(fill_color='aqua')
# 緯度/経度グリッド線を30度ごとに描画します。
map.drawmeridians(np.arange(0,360,30))
map.drawparallels(np.arange(-90,90,30))
# 定期的な緯度/経度グリッドで何らかのデータを作ります。
nlats = 73; nlons = 145; delta = 2.*np.pi/(nlons-1)
lats = (0.5*np.pi-delta*np.indices((nlats,nlons))[0,:,:])
lons = (delta*np.indices((nlats,nlons))[1,:,:])
wave = 0.75*(np.sin(2.*lats)**8*np.cos(4.*lons))
mean = 0.5*np.cos(2.*lats)*((np.sin(2.*lats))**2 + 2.)
# lat / lonグリッドのネイティブマップ投影座標を計算します。
x, y = map(lons*180./np.pi, lats*180./np.pi)
# 地図上の等高線データ。
cs = map.contour(x,y,wave+mean,15,linewidths=1.5)
plt.title('contour lines over filled continent background')
plt.show()

image.png

塗りつぶした等高線で降水量をプロット

from mpl_toolkits.basemap import Basemap, cm
# requires netcdf4-python (netcdf4-python.googlecode.com)
from netCDF4 import Dataset as NetCDFFile
import numpy as np
import matplotlib.pyplot as plt

# NWSから使用される特別な降水量カラーマップを使用して
# NWSから降雨をプロットし、ベースマップに含めます。

nc = NetCDFFile('../../../examples/nws_precip_conus_20061222.nc')
# data from http://water.weather.gov/precip/
prcpvar = nc.variables['amountofprecip']
data = 0.01*prcpvar[:]
latcorners = nc.variables['lat'][:]
loncorners = -nc.variables['lon'][:]
lon_0 = -nc.variables['true_lon'].getValue()
lat_0 = nc.variables['true_lat'].getValue()
# FigureインスタンスとAxesインスタンスを作成する
fig = plt.figure(figsize=(8,8))
ax = fig.add_axes([0.1,0.1,0.8,0.8])
# polar stereographic Basemapインスタンスを作成します。
m = Basemap(projection='stere',lon_0=lon_0,lat_0=90.,lat_ts=lat_0,\
            llcrnrlat=latcorners[0],urcrnrlat=latcorners[2],\
            llcrnrlon=loncorners[0],urcrnrlon=loncorners[2],\
            rsphere=6371200.,resolution='l',area_thresh=10000)
# 海岸線、国境、国境、地図の端を描画します。
m.drawcoastlines()
m.drawstates()
m.drawcountries()
# 緯度線を描く。
parallels = np.arange(0.,90,10.)
m.drawparallels(parallels,labels=[1,0,0,0],fontsize=10)
# 経線を描く
meridians = np.arange(180.,360.,10.)
m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10)
ny = data.shape[0]; nx = data.shape[1]
lons, lats = m.makegrid(nx, ny) # nxの緯度/経度をnxの等間隔の空間グリッドで取得します。
x, y = m(lons, lats) # map proj座標を計算します。
# 塗りつぶされた輪郭を描く
clevs = [0,1,2.5,5,7.5,10,15,20,30,40,50,70,100,150,200,250,300,400,500,600,750]
cs = m.contourf(x,y,data,clevs,cmap=cm.s3pcpn)
# add colorbar.
cbar = m.colorbar(cs,location='bottom',pad="5%")
cbar.set_label('mm')
# add title
plt.title(prcpvar.long_name+' for period ending '+prcpvar.dateofdata)
plt.show()

image.png

標高の高い高低の海面気圧の天気図をプロットする

"""
海面の圧力マップ上にHとLをプロットする
(uses scipy.ndimage.filters and netcdf4-python)
"""
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime
from mpl_toolkits.basemap import Basemap, addcyclic
from scipy.ndimage.filters import minimum_filter, maximum_filter
from netCDF4 import Dataset

def extrema(mat,mode='wrap',window=10):
    """find the indices of local extrema (min and max)
    in the input array."""
    mn = minimum_filter(mat, size=window, mode=mode)
    mx = maximum_filter(mat, size=window, mode=mode)
    # (mat == mx) true if pixel is equal to the local max
    # (mat == mn) true if pixel is equal to the local in
    # Return the indices of the maxima, minima
    return np.nonzero(mat == mn), np.nonzero(mat == mx)

# plot 00 UTC today.
date = datetime.now().strftime('%Y%m%d')+'00'

# OpenDAPデータセットを開きます。
# data=Dataset("http://nomads.ncep.noaa.gov:9090/dods/gfs/gfs/%s/gfs_%sz_anl" %\
#        (date[0:8],date[8:10]))
data=Dataset("http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd%s/gfs_hd_%sz"%\
        (date[0:8],date[8:10]))



# read lats,lons.
lats = data.variables['lat'][:]
lons1 = data.variables['lon'][:]
nlats = len(lats)
nlons = len(lons1)
# read prmsl, convert to hPa (mb).
prmsl = 0.01*data.variables['prmslmsl'][0]
# ウィンドウパラメータは、検出された高および低の数を制御する。
# (高い価値、少ない高値と低い値)
local_min, local_max = extrema(prmsl, mode='wrap', window=50)
# Basemapインスタンスを作成します。
m =\
Basemap(llcrnrlon=0,llcrnrlat=-80,urcrnrlon=360,urcrnrlat=80,projection='mill')
# 経度のラップアラウンドポイントを追加します。
prmsl, lons = addcyclic(prmsl, lons1)
# 輪郭レベル
clevs = np.arange(900,1100.,5.)
# 地図投影グリッドのx、yを見つける。
lons, lats = np.meshgrid(lons, lats)
x, y = m(lons, lats)
# 図を作成します。
fig=plt.figure(figsize=(8,4.5))
ax = fig.add_axes([0.05,0.05,0.9,0.85])
cs = m.contour(x,y,prmsl,clevs,colors='k',linewidths=1.)
m.drawcoastlines(linewidth=1.25)
m.fillcontinents(color='0.8')
m.drawparallels(np.arange(-80,81,20),labels=[1,1,0,0])
m.drawmeridians(np.arange(0,360,60),labels=[0,0,0,1])
xlows = x[local_min]; xhighs = x[local_max]
ylows = y[local_min]; yhighs = y[local_max]
lowvals = prmsl[local_min]; highvals = prmsl[local_max]
# 下の最小圧力値で、青Lとしてプロットを表示します。
xyplotted = []
# すでにdminメーター内にLまたはHがある場合はプロットしないでください。
yoffset = 0.022*(m.ymax-m.ymin)
dmin = yoffset
for x,y,p in zip(xlows, ylows, lowvals):
    if x < m.xmax and x > m.xmin and y < m.ymax and y > m.ymin:
        dist = [np.sqrt((x-x0)**2+(y-y0)**2) for x0,y0 in xyplotted]
        if not dist or min(dist) > dmin:
            plt.text(x,y,'L',fontsize=14,fontweight='bold',
                    ha='center',va='center',color='b')
            plt.text(x,y-yoffset,repr(int(p)),fontsize=9,
                    ha='center',va='top',color='b',
                    bbox = dict(boxstyle="square",ec='None',fc=(1,1,1,0.5)))
            xyplotted.append((x,y))
# 下の最大圧力値を持つ、赤色のHとしてプロットの高さ。
xyplotted = []
for x,y,p in zip(xhighs, yhighs, highvals):
    if x < m.xmax and x > m.xmin and y < m.ymax and y > m.ymin:
        dist = [np.sqrt((x-x0)**2+(y-y0)**2) for x0,y0 in xyplotted]
        if not dist or min(dist) > dmin:
            plt.text(x,y,'H',fontsize=14,fontweight='bold',
                    ha='center',va='center',color='r')
            plt.text(x,y-yoffset,repr(int(p)),fontsize=9,
                    ha='center',va='top',color='r',
                    bbox = dict(boxstyle="square",ec='None',fc=(1,1,1,0.5)))
            xyplotted.append((x,y))
plt.title('Mean Sea-Level Pressure (with Highs and Lows) %s' % date)
#          平均海面圧力(高低)
plt.show()

シェイプファイルからハリケーントラックをプロットする

"""
Cat 4または5に到達した嵐のためにAtlantic Hurricane Trackを描画します。
嵐がcat 4または5であるトラックの一部は赤く表示されます。
ESRI shapefile data from http://nationalatlas.gov/mld/huralll.html
"""
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
# ランバートコンフォーマルコニックマップ。
m = Basemap(llcrnrlon=-100.,llcrnrlat=0.,urcrnrlon=-20.,urcrnrlat=57.,
            projection='lcc',lat_1=20.,lat_2=40.,lon_0=-60.,
            resolution ='l',area_thresh=1000.)
# shapefileを読んでください。
shp_info = m.readshapefile('../../../examples/huralll020','hurrtracks',drawbounds=False)
# Cat 4に到達した嵐の名前を見つけます。
names = []
for shapedict in m.hurrtracks_info:
    cat = shapedict['CATEGORY']
    name = shapedict['NAME']
    if cat in ['H4','H5'] and name not in names:
        # 名前のついた嵐だけを使う。
        if name != 'NOT NAMED':  names.append(name)
# それらの嵐のプロットトラック。
for shapedict,shape in zip(m.hurrtracks_info,m.hurrtracks):
    name = shapedict['NAME']
    cat = shapedict['CATEGORY']
    if name in names:
        xx,yy = zip(*shape)
        # storm> Cat 4のトラックの一部を濃い赤色で表示します。
        if cat in ['H4','H5']:
            m.plot(xx,yy,linewidth=1.5,color='r')
        elif cat in ['H1','H2','H3']:
            m.plot(xx,yy,color='k')
# 海岸線、子午線および平行線を描画します。
m.drawcoastlines()
m.drawcountries()
m.drawmapboundary(fill_color='#99ffff')
m.fillcontinents(color='#cc9966',lake_color='#99ffff')
m.drawparallels(np.arange(10,70,20),labels=[1,1,0,0])
m.drawmeridians(np.arange(-100,0,20),labels=[0,0,0,1])
plt.title('Atlantic Hurricane Tracks (Storms Reaching Category 4, 1851-2004)')
plt.show()

image.png

エトポ5トポグラフィ/地形データを画像としてプロットする(指定された光源からのシェーディングの有無を問わない)。

from mpl_toolkits.basemap import Basemap, shiftgrid, cm
import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset

# etopo5地形/地形を読んでください。
url = 'http://ferret.pmel.noaa.gov/thredds/catalog/data/PMEL/etopo5.nc'
etopodata = Dataset(url)

topoin = etopodata.variables['ROSE'][:]
lons = etopodata.variables['ETOPO05_X'][:]
lats = etopodata.variables['ETOPO05_Y'][:]
# lonが20から380の代わりに-180から180になるようにデータをシフトします。
topoin,lons = shiftgrid(180.,topoin,lons,start=False)

# 地形/地形を画像としてプロットする。

# FigureインスタンスとAxesインスタンスを作成します。
fig = plt.figure()
ax = fig.add_axes([0.1,0.1,0.8,0.8])
# ベースマップの設定( 'lcc' = lambert等角錐)。
# WGS84楕円体の大小の球の半径を使用します。
m = Basemap(llcrnrlon=-145.5,llcrnrlat=1.,urcrnrlon=-2.566,urcrnrlat=46.352,\
            rsphere=(6378137.00,6356752.3142),\
            resolution='l',area_thresh=1000.,projection='lcc',\
            lat_1=50.,lon_0=-107.,ax=ax)
# nx x ny規則的に間隔を置いた5kmのネイティブ投影格子に変換する
nx = int((m.xmax-m.xmin)/5000.)+1; ny = int((m.ymax-m.ymin)/5000.)+1
topodat = m.transform_scalar(topoin,lons,lats,nx,ny)
# imshowを使って地図上にプロット画像を作成する。
im = m.imshow(topodat,cm.GMT_haxby)
# 海岸線と政治的境界を描く。
m.drawcoastlines()
m.drawcountries()
m.drawstates()
# 緯度と経度を描きます。
# 地図の左と下のラベル。
parallels = np.arange(0.,80,20.)
m.drawparallels(parallels,labels=[1,0,0,1])
meridians = np.arange(10.,360.,30.)
m.drawmeridians(meridians,labels=[1,0,0,1])
# add colorbar
cb = m.colorbar(im,"right", size="5%", pad='2%')
ax.set_title('ETOPO5 Topography - Lambert Conformal Conic')
plt.show()

# 影付きの救済プロットを作成します。

# 新しいFigure、axesインスタンスを作成します。
fig = plt.figure()
ax = fig.add_axes([0.1,0.1,0.8,0.8])
# 既存のBasemapインスタンスに新しいAxesイメージを添付します。
m.ax = ax
# 光源オブジェクトを作成します。
from matplotlib.colors import LightSource
ls = LightSource(azdeg = 90, altdeg = 20)
# 光源からの陰影を含むデータをRGB配列に変換します。
# (カラーマップを指定する必要があります)
rgb = ls.shade(topodat, cm.GMT_haxby)
im = m.imshow(rgb)
# 海岸線と政治的境界を描く。
m.drawcoastlines()
m.drawcountries()
m.drawstates()
ax.set_title('Shaded ETOPO5 Topography - Lambert Conformal Conic')
plt.show()
ARGOの位置にマーカーをプロットします。
from netCDF4 import Dataset, num2date
import time, calendar, datetime, numpy
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import urllib, os
# フォームからダウンロードしたデータ
# http://coastwatch.pfeg.noaa.gov/erddap/tabledap/apdrcArgoAll.html
filename, headers = urllib.urlretrieve('http://coastwatch.pfeg.noaa.gov/erddap/tabledap/apdrcArgoAll.nc?longitude,latitude,time&longitude>=0&longitude<=360&latitude>=-90&latitude<=90&time>=2010-01-01&time<=2010-01-08&distinct()')
dset = Dataset(filename)
lats = dset.variables['latitude'][:]
lons = dset.variables['longitude'][:]
time = dset.variables['time']
times = time[:]
t1 = times.min(); t2 = times.max()
date1 = num2date(t1, units=time.units)
date2 = num2date(t2, units=time.units)
dset.close()
os.remove(filename)
# 浮動小数点数の位置を示すマーカー付きマップ描画
m = Basemap(projection='hammer',lon_0=180)
x, y = m(lons,lats)
m.drawmapboundary(fill_color='#99ffff')
m.fillcontinents(color='#cc9966',lake_color='#99ffff')
m.scatter(x,y,3,marker='o',color='k')
plt.title('Locations of %s ARGO floats active between %s and %s' %\
        (len(lats),date1,date2),fontsize=12)
plt.show()

image.png

SSTの擬似カラープロットと海氷解析。

from mpl_toolkits.basemap import Basemap
from netCDF4 import Dataset, date2index
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime
date = datetime(2007,12,15,0) # date to plot.
# open dataset.
dataset = \
Dataset('http://www.ncdc.noaa.gov/thredds/dodsC/OISST-V2-AVHRR_agg')
timevar = dataset.variables['time']
timeindex = date2index(date,timevar) # 希望の日付の時間インデックスを見つける。
# sstを読んでください。 missing_value変数属性を使用してマスクされた配列を自動的に作成します。 
# シングルトン次元を「絞り出す」。
sst = dataset.variables['sst'][timeindex,:].squeeze()
# 氷を読む。
ice = dataset.variables['ice'][timeindex,:].squeeze()
# (グリッドボックスの中心を表す)latとlonを読む。
lats = dataset.variables['lat'][:]
lons = dataset.variables['lon'][:]
lons, lats = np.meshgrid(lons,lats)
# figure、axesインスタンスを作成します。
fig = plt.figure()
ax = fig.add_axes([0.05,0.05,0.9,0.9])
# Basemapインスタンスを作成します。
# 海岸線は使用されていないため、解像度をNoneに設定して大陸の処理をスキップします
# (これにより処理速度が向上します)
m = Basemap(projection='kav7',lon_0=0,resolution=None)
# 地図投影肢の周りに線を引く。
# 地図投影領域の色の背景。
# 土地上の欠損値がこの色を表示します。
m.drawmapboundary(fill_color='0.3')
# プロットsst、次にpcolorの氷
im1 = m.pcolormesh(lons,lats,sst,shading='flat',cmap=plt.cm.jet,latlon=True)
im2 = m.pcolormesh(lons,lats,ice,shading='flat',cmap=plt.cm.gist_gray,latlon=True)
# 平行線と子午線を描きますが、それらにラベルを付けても構いません。
m.drawparallels(np.arange(-90.,99.,30.))
m.drawmeridians(np.arange(-180.,180.,60.))
# add colorbar
cb = m.colorbar(im1,"bottom", size="5%", pad="2%")
# add a title.
ax.set_title('SST and ICE analysis for %s'%date)
plt.show()

風ベクトルと風のバーブをプロットする。

import numpy as np
import matplotlib.pyplot as plt
import datetime
from mpl_toolkits.basemap import Basemap, shiftgrid
from netCDF4 import Dataset
# specify date to plot.
yyyy=1993; mm=03; dd=14; hh=00
date = datetime.datetime(yyyy,mm,dd,hh)
# set OpenDAP server URL.
URLbase="http://nomads.ncdc.noaa.gov/thredds/dodsC/modeldata/cmd_pgbh/"
URL=URLbase+"%04i/%04i%02i/%04i%02i%02i/pgbh00.gdas.%04i%02i%02i%02i.grb2" %\
             (yyyy,yyyy,mm,yyyy,mm,dd,yyyy,mm,dd,hh)
data = Dataset(URL)
# read lats,lons
# 逆の緯度で南から北へ行く。
latitudes = data.variables['lat'][::-1]
longitudes = data.variables['lon'][:].tolist()
# 海面の圧力と10 mの風のデータを得る。
# hpの単位で0.01単位で多項式を入力します。
slpin = 0.01*data.variables['Pressure_msl'][:].squeeze()
uin = data.variables['U-component_of_wind_height_above_ground'][:].squeeze()
vin = data.variables['V-component_of_wind_height_above_ground'][:].squeeze()
# 周期的な点を手動で追加する(余弦関数を使うことができる)
slp = np.zeros((slpin.shape[0],slpin.shape[1]+1),np.float)
slp[:,0:-1] = slpin[::-1]; slp[:,-1] = slpin[::-1,0]
u = np.zeros((uin.shape[0],uin.shape[1]+1),np.float64)
u[:,0:-1] = uin[::-1]; u[:,-1] = uin[::-1,0]
v = np.zeros((vin.shape[0],vin.shape[1]+1),np.float64)
v[:,0:-1] = vin[::-1]; v[:,-1] = vin[::-1,0]
longitudes.append(360.); longitudes = np.array(longitudes)
# make 2-d grid of lons, lats
lons, lats = np.meshgrid(longitudes,latitudes)
# 正射図のベースマップを作る。
m = Basemap(resolution='c',projection='ortho',lat_0=60.,lon_0=-60.)
# 図形の作成、軸の追加
fig1 = plt.figure(figsize=(8,10))
ax = fig1.add_axes([0.1,0.1,0.8,0.8])
# 所望の輪郭レベルを設定する。
clevs = np.arange(960,1061,5)
# グリッドのネイティブx、y座標を計算します。
x, y = m(lons, lats)
# 描画する平行線と子午線を定義します。
parallels = np.arange(-80.,90,20.)
meridians = np.arange(0.,360.,20.)
# SLP輪郭をプロットする。
CS1 = m.contour(x,y,slp,clevs,linewidths=0.5,colors='k',animated=True)
CS2 = m.contourf(x,y,slp,clevs,cmap=plt.cm.RdBu_r,animated=True)
# 投影グリッド上に風ベクトルをプロットする。
# 最初にシフトグリッドを-180から180に変更します
# (経度0から360ではなく)。 さもなければ、補間は乱される。
ugrid,newlons = shiftgrid(180.,u,longitudes,start=False)
vgrid,newlons = shiftgrid(180.,v,longitudes,start=False)
# ベクトルを投影格子に変換する。
uproj,vproj,xx,yy = \
m.transform_vector(ugrid,vgrid,newlons,latitudes,31,31,returnxy=True,masked=True)
# now plot.
Q = m.quiver(xx,yy,uproj,vproj,scale=700)
# キーを作る。
qk = plt.quiverkey(Q, 0.1, 0.1, 20, '20 m/s', labelpos='W')
# 海岸線、平行線、子午線を描画します。
m.drawcoastlines(linewidth=1.5)
m.drawparallels(parallels)
m.drawmeridians(meridians)
# add colorbar
cb = m.colorbar(CS2,"bottom", size="5%", pad="2%")
cb.set_label('hPa')
# set plot title
ax.set_title('SLP and Wind Vectors '+str(date))
plt.show()

# create 2nd figure, add axes
fig2 = plt.figure(figsize=(8,10))
ax = fig2.add_axes([0.1,0.1,0.8,0.8])
# plot SLP contours
CS1 = m.contour(x,y,slp,clevs,linewidths=0.5,colors='k',animated=True)
CS2 = m.contourf(x,y,slp,clevs,cmap=plt.cm.RdBu_r,animated=True)
# マップ上に風の棘をプロットする。
barbs = m.barbs(xx,yy,uproj,vproj,length=5,barbcolor='k',flagcolor='r',linewidth=0.5)
# 海岸線、平行線、子午線を描画します。
m.drawcoastlines(linewidth=1.5)
m.drawparallels(parallels)
m.drawmeridians(meridians)
# add colorbar
cb = m.colorbar(CS2,"bottom", size="5%", pad="2%")
cb.set_label('hPa')
# set plot title.
ax.set_title('SLP and Wind Barbs '+str(date))
plt.show()

NYとロンドンの間に大きな円(大圏コース)を描く。

from mpl_toolkits.basemap import Basemap
import numpy as np
import matplotlib.pyplot as plt
# 新しいFigure、Axesのインスタンスを作成します。
fig=plt.figure()
ax=fig.add_axes([0.1,0.1,0.8,0.8])
# セットアップメルケータマップ投影。
m = Basemap(llcrnrlon=-100.,llcrnrlat=20.,urcrnrlon=20.,urcrnrlat=60.,\
            rsphere=(6378137.00,6356752.3142),\
            resolution='l',projection='merc',\
            lat_0=40.,lon_0=-20.,lat_ts=20.)
# nylat, nylon are lat/lon of New York
nylat = 40.78; nylon = -73.98
# lonlat, lonlon are lat/lon of London.
lonlat = 51.53; lonlon = 0.08
# NYとロンドンの間に大きな円(大圏コース)ルートを描く
m.drawgreatcircle(nylon,nylat,lonlon,lonlat,linewidth=2,color='b')
m.drawcoastlines()
m.fillcontinents()
# draw parallels
m.drawparallels(np.arange(10,90,20),labels=[1,1,0,1])
# draw meridians
m.drawmeridians(np.arange(-180,180,30),labels=[1,1,0,1])
ax.set_title('Great Circle from New York to London')
plt.show()

image.png

昼夜のターミネーターを地図上に描く。

import numpy as np
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from datetime import datetime
# ミラー投影
map = Basemap(projection='mill',lon_0=180)
# 海岸をプロットし、ラベルの子午線と平行を描きます。
map.drawcoastlines()
map.drawparallels(np.arange(-90,90,30),labels=[1,0,0,0])
map.drawmeridians(np.arange(map.lonmin,map.lonmax+30,60),labels=[0,0,0,1])
# 大陸の 'サンゴ'(zorder = 0)、カラーウェットエリア 'アクア'
map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='coral',lake_color='aqua')
# 夜間領域を陰影付けし、地図が表示されるようにアルファ透明度を設定します。 
# UTCで現在の時刻を使用します。
date = datetime.utcnow()
CS=map.nightshade(date)
plt.title('Day/Night Map for %s (UTC)' % date.strftime("%d %b %Y %H:%M:%S"))
plt.show()

image.png

15
21
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
15
21

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?