2017-10-02 6 views
1

緯度経度が-180から+180にスワップするような、緯度/経度/データのペアがアンティメリディアを横切る場合、どのようにしてpcolor(mesh)地球全体を満たすグリッドセルを描くことから?私の問題は、basemapではなくcartopyを使用していることを除いて、hereと同じです。リンクされた質問(約basemap)に近い5年前のコメントは、cartopy解決策があると主張しているが、これは掲載されていない。グリッドのないグリッド線(水平線)を防止する

例コード:

#!/usr/bin/env python3.6 

import numpy 
import matplotlib.pyplot 
import cartopy.crs 

lons = numpy.array([[-174.719, -175.297, -175.883], 
     [-175.164, -175.734, -176.312], 
     [-175.594, -176.164, -176.734], 
     [-176.016, -176.578, -177.148], 
     [-176.43 , -176.984, -177.547], 
     [-176.836, -177.383, -177.938], 
     [-177.227, -177.773, -178.312], 
     [-177.609, -178.148, -178.688], 
     [-177.984, -178.516, -179.047], 
     [-178.352, -178.875, -179.398], 
     [-179.727, 179.766, 179.266], 
     [ 179.945, 179.445, 178.945], 
     [ 179.625, 179.133, 178.641], 
     [ 179.312, 178.828, 178.336], 
     [ 179.008, 178.523, 178.039], 
     [ 178.711, 178.234, 177.75 ], 
     [ 178.414, 177.945, 177.469], 
     [ 178.133, 177.656, 177.188], 
     [ 177.844, 177.383, 176.914], 
     [ 177.57 , 177.109, 176.648]]) 

lats = numpy.array([[ 67.391, 67.492, 67.586], 
     [ 67.055, 67.148, 67.25 ], 
     [ 66.711, 66.812, 66.906], 
     [ 66.375, 66.469, 66.562], 
     [ 66.031, 66.125, 66.219], 
     [ 65.688, 65.781, 65.875], 
     [ 65.344, 65.438, 65.523], 
     [ 65. , 65.094, 65.18 ], 
     [ 64.656, 64.742, 64.836], 
     [ 64.312, 64.398, 64.484], 
     [ 62.922, 63. , 63.086], 
     [ 62.57 , 62.648, 62.734], 
     [ 62.219, 62.297, 62.383], 
     [ 61.867, 61.945, 62.023], 
     [ 61.516, 61.594, 61.672], 
     [ 61.164, 61.242, 61.32 ], 
     [ 60.812, 60.891, 60.961], 
     [ 60.812, 60.891, 60.961], 
     [ 60.461, 60.531, 60.609], 
     [ 60.102, 60.18 , 60.25 ]]) 

data = numpy.array([[ 231.73, 231.56, 231.22], 
     [ 231.72, 231.72, 231.72], 
     [ 232.24, 232.73, 233.37], 
     [ 233.22, 233.69, 234.01], 
     [ 234.33, 234.94, 235.39], 
     [ 234.5 , 235.11, 235.71], 
     [ 235.41, 235.71, 236. ], 
     [ 235.27, 235.72, 236.31], 
     [ 234.67, 235.43, 235.73], 
     [ 235.43, 236.17, 235.88], 
     [ 236.18, 236.18, 236.18], 
     [ 236.07, 236.36, 236.79], 
     [ 235.8 , 236.1 , 235.8 ], 
     [ 236.84, 236.84, 236.55], 
     [ 238.27, 238.27, 238.54], 
     [ 237.72, 237.44, 237.72], 
     [ 238.42, 238.28, 238.28], 
     [ 238.57, 238.57, 238.43], 
     [ 240.17, 240.04, 239.65], 
     [ 241.21, 241.21, 241.09]]) 

proj = cartopy.crs.Mollweide() 
ax = matplotlib.pyplot.axes(projection=proj) 
trans = proj.transform_points(cartopy.crs.Geodetic(), lons, lats) 
ax.coastlines() 
ax.pcolormesh(trans[:, :, 0], trans[:, :, 1], data, transform=proj) 

matplotlib.pyplot.savefig("/tmp/test.png") 

の予想される出力は、北太平洋のどこかに中心にデータのビットを持つマップになります。現実には、私は、地球全体の幅にまたがる非常に細長いマップを取得:

Elongated map

私はより簡単に質問にそれを組み込むことができるようなポイントの数が少ないにデータを制限しました、実際には私は両極を常に横切る極衛星データの完全軌道を持っているので、常にアンチメリディアンを横切っています。実際の軌道のための結果は次のようになります。中央経度を変更

Map full of nonsense

は、問題を再配置します。私は地図の端を横切る場所から中心の経度を選ぶことによって重大度を減らすことができます。この例では、前のマップと同じデータがプロットされているが、90°Eの中心経度と:

Still bad but different

This pull request 2012から関連する表示され、そう明らかに、関連する機能があることが想定されますしかし、私はそれをどのように使用するのか手がかりがありません。この問題は、どのグローバルマップ投影でも発生します。私はcartopy 0.15.1を使用しています。

これを正しくプロットするにはどうすればよいですか?

+0

['Mollweide'投影](http://scitools.org.uk/cartopy/docs/latest/crs/projections.html#mollweide)にはオプション' central_longitude'があります。あなたはその価値を変えようとしましたか? – tom

+0

@tom良い質問です。それは解決しませんが、手がかりを提供します。これは、緯度/経度座標の折り返しではなく、投影された座標の折り返しです。更新された地図を見る。 – gerrit

答えて

0

編集:この上の2つのフォローアップの質問があることを注:*総ショーで、与えられた問題の回避策を希望

*


元の答え:

私は次のことが大きな助けになるかどうかはわかりませんが、試してみましょう。 -180と180の間のシフトが発生する位置でデータを分割することができれば、2つの異なるpcolorプロットをプロットすることができ、このようにしてpcolorshapesが地球周りを1周回するのを防ぎます。

テストデータを使用すると、これは簡単です。私はちょうどデータポイントの1つの記号を変えました。しかし、1行分のデータを残すこともできます。次にデータをプロットすることで、所望の結果が得られます。両方のプロットでカラーマップが同期していることを確認するためには、正規化が必要です。

norm = plt.Normalize(data.min(), data.max()) 
    ax.pcolormesh(trans[:10, :, 0], trans[:10, :, 1], data[:10,:], transform=proj, norm=norm) 
    ax.pcolormesh(trans[10:, :, 0], trans[10:, :, 1], data[10:,:], transform=proj, norm=norm) 

enter image description here ..... enter image description here

完全なコード:

まず
import numpy 
import matplotlib.pyplot as plt 
import cartopy.crs 

lons = numpy.array([[-174.719, -175.297, -175.883], 
     [-175.164, -175.734, -176.312], 
     [-175.594, -176.164, -176.734], 
     [-176.016, -176.578, -177.148], 
     [-176.43 , -176.984, -177.547], 
     [-176.836, -177.383, -177.938], 
     [-177.227, -177.773, -178.312], 
     [-177.609, -178.148, -178.688], 
     [-177.984, -178.516, -179.047], 
     [-178.352, -178.875, -179.398], 
     [ 179.999, 179.766, 179.266], #<- changed sign here 
     [ 179.945, 179.445, 178.945], 
     [ 179.625, 179.133, 178.641], 
     [ 179.312, 178.828, 178.336], 
     [ 179.008, 178.523, 178.039], 
     [ 178.711, 178.234, 177.75 ], 
     [ 178.414, 177.945, 177.469], 
     [ 178.133, 177.656, 177.188], 
     [ 177.844, 177.383, 176.914], 
     [ 177.57 , 177.109, 176.648]]) 

lats = numpy.array([[ 67.391, 67.492, 67.586], 
     [ 67.055, 67.148, 67.25 ], 
     [ 66.711, 66.812, 66.906], 
     [ 66.375, 66.469, 66.562], 
     [ 66.031, 66.125, 66.219], 
     [ 65.688, 65.781, 65.875], 
     [ 65.344, 65.438, 65.523], 
     [ 65. , 65.094, 65.18 ], 
     [ 64.656, 64.742, 64.836], 
     [ 64.312, 64.398, 64.484], 
     [ 62.922, 63. , 63.086], 
     [ 62.57 , 62.648, 62.734], 
     [ 62.219, 62.297, 62.383], 
     [ 61.867, 61.945, 62.023], 
     [ 61.516, 61.594, 61.672], 
     [ 61.164, 61.242, 61.32 ], 
     [ 60.812, 60.891, 60.961], 
     [ 60.812, 60.891, 60.961], 
     [ 60.461, 60.531, 60.609], 
     [ 60.102, 60.18 , 60.25 ]]) 

data = numpy.array([[ 231.73, 231.56, 231.22], 
     [ 231.72, 231.72, 231.72], 
     [ 232.24, 232.73, 233.37], 
     [ 233.22, 233.69, 234.01], 
     [ 234.33, 234.94, 235.39], 
     [ 234.5 , 235.11, 235.71], 
     [ 235.41, 235.71, 236. ], 
     [ 235.27, 235.72, 236.31], 
     [ 234.67, 235.43, 235.73], 
     [ 235.43, 236.17, 235.88], 
     [ 236.18, 236.18, 236.18], 
     [ 236.07, 236.36, 236.79], 
     [ 235.8 , 236.1 , 235.8 ], 
     [ 236.84, 236.84, 236.55], 
     [ 238.27, 238.27, 238.54], 
     [ 237.72, 237.44, 237.72], 
     [ 238.42, 238.28, 238.28], 
     [ 238.57, 238.57, 238.43], 
     [ 240.17, 240.04, 239.65], 
     [ 241.21, 241.21, 241.09]]) 
print lons.shape, lats.shape, data.shape 
proj = cartopy.crs.Mollweide() 
ax = plt.axes(projection=proj) 
trans = proj.transform_points(cartopy.crs.Geodetic(), lons, lats) 
ax.coastlines() 
norm = plt.Normalize(data.min(), data.max()) 
ax.pcolormesh(trans[:10, :, 0], trans[:10, :, 1], data[:10,:], transform=proj, norm=norm) 
ax.pcolormesh(trans[10:, :, 0], trans[10:, :, 1], data[10:,:], transform=proj, norm=norm) 

plt.show() 
+0

これはテストデータには有効ですが、スキャンラインがアンチメリディアンと交差すると失敗します。つまり、[-179.8,179.1,179.3]です。回避策はおそらく構築できると思うが、カートリスト内にネイティブな解決策があるのだろうかと思う。 – gerrit

+0

マスクされた配列を使用することができます(最初のプロットではすべて負の値、2番目以降の場合はすべて正の値をマスクします)。問題は、(非単調なグリッドで使用するように設計されていない)matplotlibのpcolormeshの関数に与えられたデータによって引き起こされるように私は強く疑いがcartopy内からの溶液です。 – ImportanceOfBeingErnest

+0

リンク先の質問にpelsonがあることが示され、「pcolormeshのラップアラウンドを処理する」というプルリクエストがあるので、私はネイティブのcartopyソリューションを望んでいました。私はマスキングの方法を試しましたが、何らかの理由で、 'ma.masked_where(lons> 0、bts)'をプロットするときにストライプを取得していて、私が間違っていることを理解しようとしています。 – gerrit

1

、再現するためのいくつかのデータとコードの一部を提供するためのおかげで - それは私がすぐにできたことを意味問題自体を重視し、問題を再現することはできません。

カートリストとベースマップの主な違いは、カートリストがベクター/ラスタ変換を処理できることです。カートリッジをベースマップの方法で動作させることは、ユーザーがデータを自分自身で変換することが重要です。あなたが提供した例は、lats/lonsを手動でターゲット投影に変換することによってこれを正確に行います。多大な注意を払わなければ、遭遇したようなアンチメリディアンの問題をすばやく見つけることになります。ありがたいことに、cartopy にデータ変換に関して非常に注意を払っており、私はあなたにそれを利用することをお勧めします。擬似コードで

、あなたのコードはありません:

create a mollweide map 
convert your lats/lons to mollweide coordinate system 
plot newly converted mollweide data on mollweide map 

は、実際には、我々はcartopyとのパラダイムを変更すると実行します。

create a mollweide map 
plot lat/lon data on mollweide map 

そうすることによって、我々は必要なcartopy与えていますコンテキストを使用してデータを正しく変換します。私はPlateCarree投影を使用していませんでした。この例では

ax.pcolormesh(lons, lats, data, transform=ccrs.PlateCarree()) 

あなたのコードに大きな変更は(ラッツ/ロンでの)元のデータ、あなたが手で変換されないの座標をプロットすることです我々は現在、測地的なpcolormeshボックス(つまり大きな円)を実装しておらず、本質的に一定の緯度/経度のボックスを生成しているので、測地座標系です。

これを使用すると、あなたの質問の最初の画像に非常によく似たプロットが作成されます。その理由は、定義しているボックスの中にはPlateCarreeの投影されたスペースの幅が〜360度の(平らな紙ですが、ラップアラウンド/アンチメリディアンは何も知りません)という幅があります。

ここでは、意図した例を見てみましょう。あなたは測地線の観点から考えるならば、あなたは次のコードは、マップのいずれかの側にある2つの小箱を生産するために期待するかもしれない:

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

box = sgeom.box(minx=170, maxx=-170, miny=40, maxy=60) 

proj = ccrs.Mollweide() 

ax = plt.axes(projection=proj) 
ax.coastlines() 
ax.add_geometries([box], ccrs.PlateCarree(), facecolor='coral', 
        edgecolor='black', alpha=0.5) 

plt.show() 

Big box, when we wanted two little ones...

悲しいかなが、それは我々が得るものではありません。これは、Plate Carree投影が2点のデカルト投影であり、2点間の唯一の有効な線が直線であることを覚えていれば理にかなっています - それはアンチメリディアンを包むことについて何も知らない。

は(それは注目に値する:私たちは、その後、我々は与えられた点の間で大きな円を描く測地ジオメトリ投影を変更し、所望のボックスを取得した場合)

だから我々は、ボックスの座標を必要とする目的のボックスを生成するために360度に近づくものではなく、小さなx範囲を持つこと。ありがたいことにcartopyは、180度を越えるPlateCarree座標値を定義することを可能にします - これは小さなx範囲のPlateCarreeボックスを定義するための鍵です。

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

box = sgeom.box(minx=170, maxx=190, miny=40, maxy=60) 

proj = ccrs.Mollweide() 

ax = plt.axes(projection=proj) 
ax.coastlines() 
ax.add_geometries([box], ccrs.PlateCarree(), facecolor='coral', 
        edgecolor='black', alpha=0.5) 

two boxes - hurrah!

だから、戻ってあなたの例に行く - 私たちは本当に測地パッチを定義している/緯度ロン、の束を持っています。 Cartopyはまだ測地座標をpcolormeshできません - 回避策はpcolormeshになりますPlateCarree座標。測地座標とPlateCarree座標ポイントは互換性がありますが、それらは基本的に異なるトポロジーを持っています。

この例では、0より下の値に360を追加することで、データを有効なPlateCarreeトポロジに変換することができます。残念ながら、これは中心子午線を横切るジオメトリでは機能しません。関与しており、IMOをカートリピアに拡張するのに便利です。

最終的なコードは、今のようになります。興味のある方

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

lons = np.array([[-174.719, -175.297, -175.883], 
     [-175.164, -175.734, -176.312], 
     [-175.594, -176.164, -176.734], 
     [-176.016, -176.578, -177.148], 
     [-176.43 , -176.984, -177.547], 
     [-176.836, -177.383, -177.938], 
     [-177.227, -177.773, -178.312], 
     [-177.609, -178.148, -178.688], 
     [-177.984, -178.516, -179.047], 
     [-178.352, -178.875, -179.398], 
     [-179.727, 179.766, 179.266], 
     [ 179.945, 179.445, 178.945], 
     [ 179.625, 179.133, 178.641], 
     [ 179.312, 178.828, 178.336], 
     [ 179.008, 178.523, 178.039], 
     [ 178.711, 178.234, 177.75 ], 
     [ 178.414, 177.945, 177.469], 
     [ 178.133, 177.656, 177.188], 
     [ 177.844, 177.383, 176.914], 
     [ 177.57 , 177.109, 176.648]]) 

lats = np.array([[ 67.391, 67.492, 67.586], 
     [ 67.055, 67.148, 67.25 ], 
     [ 66.711, 66.812, 66.906], 
     [ 66.375, 66.469, 66.562], 
     [ 66.031, 66.125, 66.219], 
     [ 65.688, 65.781, 65.875], 
     [ 65.344, 65.438, 65.523], 
     [ 65. , 65.094, 65.18 ], 
     [ 64.656, 64.742, 64.836], 
     [ 64.312, 64.398, 64.484], 
     [ 62.922, 63. , 63.086], 
     [ 62.57 , 62.648, 62.734], 
     [ 62.219, 62.297, 62.383], 
     [ 61.867, 61.945, 62.023], 
     [ 61.516, 61.594, 61.672], 
     [ 61.164, 61.242, 61.32 ], 
     [ 60.812, 60.891, 60.961], 
     [ 60.812, 60.891, 60.961], 
     [ 60.461, 60.531, 60.609], 
     [ 60.102, 60.18 , 60.25 ]]) 

data = np.array([[ 231.73, 231.56, 231.22], 
     [ 231.72, 231.72, 231.72], 
     [ 232.24, 232.73, 233.37], 
     [ 233.22, 233.69, 234.01], 
     [ 234.33, 234.94, 235.39], 
     [ 234.5 , 235.11, 235.71], 
     [ 235.41, 235.71, 236. ], 
     [ 235.27, 235.72, 236.31], 
     [ 234.67, 235.43, 235.73], 
     [ 235.43, 236.17, 235.88], 
     [ 236.18, 236.18, 236.18], 
     [ 236.07, 236.36, 236.79], 
     [ 235.8 , 236.1 , 235.8 ], 
     [ 236.84, 236.84, 236.55], 
     [ 238.27, 238.27, 238.54], 
     [ 237.72, 237.44, 237.72], 
     [ 238.42, 238.28, 238.28], 
     [ 238.57, 238.57, 238.43], 
     [ 240.17, 240.04, 239.65], 
     [ 241.21, 241.21, 241.09]]) 

proj = ccrs.PlateCarree(central_longitude=180) 
ax = plt.axes(projection=proj) 
ax.coastlines('50m') 
ax.margins(0.3) 

lons[lons < 0] += 360 
ax.pcolormesh(lons, lats, data, transform=ccrs.PlateCarree()) 

plt.show() 

the resulting geometry

、私は一般的に、プレートカレのものに測地pcolormeshの境界を変換する関数を追加するcartopy機能要求を開くことをお勧めします。 cartopyトラッカーはhttps://github.com/SciTools/cartopy/issues/newにあります。

関連する問題