2016-04-04 9 views
4

いくつかの衛星ベースの地球観測製品は緯度/経度情報を提供し、他のものは所定のグリッド投影内にX/Y座標を提供します(例を参照)。 2番目の場合の私のアプローチは、与えられたX/Y値がベースマップの座標と等しくなるようにデータプロバイダによって与えられたものと同じパラメータ(射影、楕円体、マップの原点)を持つBasemapマップを設定することです。しかし、私がそうするなら、地理位置はBasemapの海岸線を含む他のデータセットと一致しません。 私はこれを、異なる信頼できるソースからの3つの異なるデータセットで経験しました。最小の例として、私はSouth polar StereographicグリッドのX/Y座標とイメージの四隅すべての対応する緯度/経度座標の両方を含むU.S. Geological Surveyによって提供されるLandsatデータを使用します。Basemap南極立体地図プロジェクション座標が同じ投影のデータセットの座標と一致しないのはなぜですか?

我々は(ID:LC82171052016079LGN00):取得ランドサットメタファイルから

CORNER_UL_LAT_PRODUCT = -66.61490 CORNER_UL_LON_PRODUCT = -61.31816 CORNER_UR_LAT_PRODUCT = -68.74325 CORNER_UR_LON_PRODUCT = -58.04533 CORNER_LL_LAT_PRODUCT = -67.68721 CORNER_LL_LON_PRODUCT = -67.01109 CORNER_LR_LAT_PRODUCT = -69.94052 CORNER_LR_LON_PRODUCT = -64.18581 CORNER_UL_PROJECTION_X_PRODUCT = -2259300.000 CORNER_UL_PROJECTION_Y_PRODUCT = 1236000.000 CORNER_UR_PROJECTION_X_PRODUCT = -1981500.000 コーナー_UR_PROJECTION_Y_PRODUCT = 1236000.000 CORNER_LL_PROJECTION_X_PRODUCT = -2259300.000 CORNER_LL_PROJECTION_Y_PRODUCT = 958500.000 CORNER_LR_PROJECTION_X_PRODUCT = -1981500.000 CORNER_LR_PROJECTION_Y_PRODUCT = 958500.000

...

GROUP = PROJECTION_PARAMETERS MAP_PROJECTION = "PS" DATUM = "WGS84" FALSE_EASTING = 0 FALSE_NORTHING = 0 GRID_CELL_SIZE_PANCHROMATIC = 15.00 GRI D_CELL_SIZE_REFLECTIVE = 30.00 GRID_CELL_SIZE_THERMAL = 30.00 ORIENTATION = "NORTH_UP" RESAMPLING_OPTION = "CUBIC_CONVOLUTION" END_GROUP = PROJECTION_PARAMETERS

右の地図投影でベースマップを使用することにより、我々は緯度/経度値からコーナーを導出することができるはずですX/Y値:

import numpy as np 
from mpl_toolkits.basemap import Basemap 

m=Basemap(resolution='h',projection='spstere', ellps='WGS84', boundinglat=-60,lon_0=180, lat_ts=-71) 

x_crn=np.array([-2259300,-1981500,-2259300,-1981500])# upper left, upper right, lower left, lower right 
y_crn=np.array([1236000, 1236000, 958500, 958500])# upper left, upper right, lower left, lower right 

x0, y0= m(0, -90) 
#Basemap coordinates at the south pole 
#note that (0,0) of the Basemap is in a corner of the map, 
#while other data sets use the south pole. 
#This is easy to take into account: 
lon_crn, lat_crn = m(x0-x_crn, y0-y_crn, inverse=True) 


print 'lon_crn: '+str(lon_crn) 
print 'lat_crn: '+str(lat_crn) 

返します

lon_crn: [-61.31816102 -58.04532791 -67.01108782 -64.1858106 ] 
lat_crn: [-67.23548626 -69.3099076 -68.28071626 -70.47651326] 

見ることができるように、経度は、メタファイルからのもので特定の精度に一致しますが、緯度は低くなります。

私はで緯度を近似することができる。

lat_crn=(lat_crn+90.)*1.0275-90. 

しかし、これは本当に満足ではありません。

これはX/Yのコーナーを使用するメタファイル(赤ベースマップdrawcoastlines())から座標場合に画像が配置されている方法である。 enter image description here 、これはそれが緯度/経度角を使用してどのように見えるかである: enter image description here

この場合、私は単純にlat/lon座標を使うことができますが、前述のようにX/Y座標のみで提供されるデータセット(thisなど)があります。これはBasemap投影に頼ることを非常に重要にします。私は潜在的な回避策としてデータを再投影するための他のモジュールがあることを知っていますが、他のモジュールなしで動作し、再投影によってエラーが発生する可能性があります。

この問題はさまざまなデータセットで発生するため、私はBasemapモジュールのバグだと思っていますが、何度も何度も同じミスを繰り返したり間違った予想をしている可能性があります。

+0

mpl_toolkits.basemapのバージョンは1.0.7 – Dusch

+0

楕円体に問題がありますか?私は間違ったことは見ませんが、緯度が違うという事実はこれを示しています... – sulkeh

+0

これは正しい@sulkehですが、Basemapの初期化のために他の楕円体をテストすると数キロメートルの北/我々が直面している〜60kmシフトへ)。 WGS84と比較して、CMPというもの(Comm。des Poids et Mesures 1799)とplessis(Plessis 1817)をパラメータとして使ってみました。しかし、楕円体のベースマップ処理には共通のエラーが存在する可能性があります。 – Dusch

答えて

2

私はいくつかの実験を行いました。変化するようです。lat_tsは、projection='spstere'で効果がありません。実際には、プロジェクションの緯度は、割り当てられる値に関係なく暗黙的にlat_ts=-90.と仮定されているようです。

次のように、あなたの例では、ベースマップを構築なるように、私は、代わりにprojection='stere'を使用してより多くの成功を持っていた:

m=Basemap(width=5400000., height=5400000., projection='stere', 
      ellps='WGS84', lon_0=180., lat_0=-90., lat_ts=-71.) 

あなたは緯度と経度のコーナーの代わりに、幅と高さを設定することを好むことあなたのアプリケーションのためのプロットの。

関連する問題