15

私は次の電子メールを受け取り、この質問への回答がすべての人に利用可能であることを確認したいと思いました。


やあ、

日付線を横切り、左側に東アジア、右側に北アメリカの西を示すカートピーを使用して、単純な緯度経度マップを設定したいと思います。次のグーグルマップは大まかに私が求めているものです:

https://maps.google.co.uk/?ll=56.559482,-175.253906&spn=47.333523,133.066406&t=m&z=4

グーグルマップのスクリーンショット

これはCartopyで実行できますか?

4

2 に答える 2

37

良い質問。これはおそらく何度も出てくるものなので、実際にあなたの特定の質問に答える前に、このステップバイステップを見ていきます。今後の参考のために、次の例はcartopyv0.5で作成されました。

まず、デフォルトの「緯度経度」(またはより技術的にはPlateCarree)投影は、-180〜180の前方範囲で機能することに注意することが重要です。これは、標準のPlateCarree投影をこれを超えてプロットできないことを意味します。これにはいくつかの理由がありますが、そのほとんどは、ベクトルとラスターの両方を投影するときにカートピーがより多くの作業を行う必要があるという事実に要約されますたとえば、単純な海岸線)。残念ながら、作成しようとしているプロットには、まさにこの機能が必要です。この制限を写真に入れるために、デフォルトのPlateCarreeプロジェクションは次のようになります。

import cartopy.crs as ccrs
import matplotlib.pyplot as plt

proj = ccrs.PlateCarree(central_longitude=0)

ax1 = plt.axes(projection=proj)
ax1.stock_img()
plt.title('Global')

plt.show()

グローバルプレートキャリー

このマップに描画できる単一の長方形は、合法的に拡大された領域にすることができます(ここには少し高度なコードがありますが、画像は1000語の価値があります)。

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import shapely.geometry as sgeom

box = sgeom.box(minx=-90, maxx=45, miny=15, maxy=70)
x0, y0, x1, y1 = box.bounds

proj = ccrs.PlateCarree(central_longitude=0)

ax1 = plt.subplot(211, projection=proj)
ax1.stock_img()
ax1.add_geometries([box], proj, facecolor='coral', 
                   edgecolor='black', alpha=0.5)
plt.title('Global')

ax2 = plt.subplot(212, projection=proj)
ax2.stock_img()
ax2.set_extent([x0, x1, y0, y1], proj)
plt.title('Zoomed in area')

plt.show()

2番目に拡大されたマップでグローバル

残念ながら、必要なプロットには、この投影法で2つの長方形が必要になります。

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import shapely.geometry as sgeom

box = sgeom.box(minx=120, maxx=260, miny=15, maxy=80)

proj = ccrs.PlateCarree(central_longitude=0)

ax1 = plt.axes(projection=proj)
ax1.stock_img()
ax1.add_geometries([box], proj, facecolor='coral', 
                   edgecolor='black', alpha=0.5)
plt.title('Target area')

plt.show()

グローバルプロット上のターゲット長方形

したがって、標準のPlateCarree定義を使用して、日付線と交差するマップを描画することはできません。代わりに、PlateCarree定義の中心経度を変更して、ターゲットとする領域の単一のボックスを描画できるようにすることができます。

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import shapely.geometry as sgeom

box = sgeom.box(minx=120, maxx=260, miny=15, maxy=80)
x0, y0, x1, y1 = box.bounds

proj = ccrs.PlateCarree(central_longitude=180)
box_proj = ccrs.PlateCarree(central_longitude=0)

ax1 = plt.subplot(211, projection=proj)
ax1.stock_img()
ax1.add_geometries([box], box_proj, facecolor='coral', 
                   edgecolor='black', alpha=0.5)
plt.title('Global')

ax2 = plt.subplot(212, projection=proj)
ax2.stock_img()
ax2.set_extent([x0, x1, y0, y1], box_proj)
plt.title('Zoomed in area')

plt.show()

中心経度が180の2つのPlateCarreプロット、1つはグローバル、もう1つはターゲット領域にズームイン

うまくいけば、それはあなたがあなたの目標マップを達成するためにあなたがしなければならないことをあなたに示すでしょ、上のコードはあなたの目標を達成するために少し複雑かもしれません、それで少し単純化するために、あなたが望むプロットを生成するために私が書くコードは次のようになります:

import cartopy.feature
import cartopy.crs as ccrs
import matplotlib.pyplot as plt    

ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
ax.set_extent([120, 260, 15, 80], crs=ccrs.PlateCarree())

# add some features to make the map a little more polished
ax.add_feature(cartopy.feature.LAND)
ax.add_feature(cartopy.feature.OCEAN)
ax.coastlines('50m')

plt.show()

最終的なターゲットマップ

これは長い答えでした。うまくいけば、私は質問に答えただけでなく、マップ作成とカートピーのより複雑な詳細のいくつかをより明確にして、将来の問題をスムーズにするのに役立てることができます。

乾杯、

于 2012-12-13T09:58:34.177 に答える
0

上記のベンジミンのコメントの詳細については、

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from matplotlib.ticker import AutoMinorLocator, FixedLocator, MultipleLocator

def map_common(ax1,gl_loc=[True,True,False,True],gl_lon_info=range(-180,180,60),gl_dlat=30):

    ax1.coastlines(color='silver',linewidth=1.)

    gl = ax1.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                       linewidth=0.6, color='gray', alpha=0.5, linestyle='--')

    gl.ylabels_left = gl_loc[0]
    gl.ylabels_right = gl_loc[1]
    gl.xlabels_top = gl_loc[2]
    gl.xlabels_bottom = gl_loc[3]

    gl.xlocator = FixedLocator(gl_lon_info)
    gl.ylocator = MultipleLocator(gl_dlat)
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER
    gl.xlabel_style = {'size': 11, 'color': 'k'}
    gl.ylabel_style = {'size': 11, 'color': 'k'}


lon_boundary=np.arange(-240,-60,1.)
lat_boundary=np.arange(15,75,1.)
data=np.ones([lat_boundary.shape[0]-1,lon_boundary.shape[0]-1]) ## Data dimension is 1 less than boundaries
data=data*lat_boundary[:-1,None]

lon_offset=-150  ##
x,y=np.meshgrid(lon_boundary-lon_offset,lat_boundary)

fig=plt.figure()
fig.set_size_inches(7.5,5) ## (xsize, ysize)
ax1=fig.add_subplot(111,projection=ccrs.PlateCarree(central_longitude=lon_offset))
ax1.set_extent([-250,-50,10,80],crs=ccrs.PlateCarree())

props=dict(vmin=0,vmax=90,cmap=plt.cm.get_cmap('bone'),alpha=0.8)
cs=ax1.pcolormesh(x,y,data,**props)
ax1.set_title('Lon_Offset=-90')
map_common(ax1,gl_lon_info=[-180,-120,-60,120,],gl_dlat=15)

fnout='./map_over_dateline.png'
#plt.show()
plt.savefig(fnout,bbox_inches='tight',dpi=150)

このプログラムの出力

于 2019-12-05T16:25:03.040 に答える