2016-10-08 3 views
1

matplotlibとcartopyを使用してChoroplethマップを作成しようとしていますが、シェープファイルを最初にプロットする必要があります。しかし、同様の質問にherehereと尋ねられても、私はそうすることはできませんでした。私は投影または境界が誤って指定されていると思われます。カートファイルの形状ファイルの一致投影

マイシェープ投影

PROJCS["WGS_1984_UTM_Zone_32Nz", 
    GEOGCS["GCS_WGS_1984", 
     DATUM["WGS_1984", 
      SPHEROID["WGS_84",6378137,298.257223563]], 
     PRIMEM["Greenwich",0], 
     UNIT["Degree",0.017453292519943295]], 
    PROJECTION["Transverse_Mercator"], 
    PARAMETER["False_Easting",32500000], 
    PARAMETER["False_Northing",0], 
    PARAMETER["Central_Meridian",9], 
    PARAMETER["Scale_Factor",0.9996], 
    PARAMETER["Latitude_Of_Origin",0], 
    UNIT["Meter",1]] 

を有し、私が話している場合、約vg250_2010-01-01.utm32w.shape.ebenen/vg250_ebenen-historisch/de1001/vg250_gem.shp

私のコードは、私は、対応する方法を使用して境界を得

#!/usr/local/bin/python 
# -*- coding: utf8 -*- 

import cartopy.crs as ccrs 
import cartopy.io.shapereader as shpreader 
import matplotlib.pyplot as plt 

fname = 'path/vg250_gem.shp' 
proj = ccrs.TransverseMercator(central_longitude=0.0,central_latitude=0.0, 
          false_easting=32500000.0,false_northing=0.0, 
          scale_factor=0.9996) 
municipalities = list(shpreader.Reader(fname).geometries()) 
ax = plt.axes(projection=proj) 
plt.title('Deutschland') 
ax.add_geometries(municipalities,proj,edgecolor='black',facecolor='gray',alpha=0.5) 
ax.set_extent([32458044.649189778*0.9, 5556418.748046352*1.1, 32465287.307457082*0.9, 5564153.5456742775*1.1],proj) 
plt.show() 

あるhereをダウンロードすることができフィオナから。 Pythonがエラーをスローする

Traceback (most recent call last): 
    File "***/src/analysis/test.py", line 16, in <module> 
ax.set_extent([32458044.649189778, 5556418.748046352, 32465287.307457082, 5564153.5456742775],proj) 
    File "/usr/local/lib/python2.7/site-packages/cartopy/mpl/geoaxes.py", line 652, in set_extent 
ylim=self.projection.y_limits)) 
ValueError: Failed to determine the required bounds in projection coordinates. Check that the values provided are within the valid range (x_limits=[-20000000.0, 20000000.0], y_limits=[-10000000.0, 10000000.0]). 
[Finished in 53.9s with exit code 1] 

これは私には意味がありません。また、ccrs.UTM()を使って実験すると、白い部分を示すプロットが得られます。もし誰かが私にこの問題を解決する方法を教えてもらえると感謝します。ありがとうございました!

+0

私はcartopyなどについて何も知りませんが、エラーメッセージに次のツー最後の行を読んで、いくつかの必要なxと記載されているyの限界が、set_extent 'へのお電話での限界があります'(エラーの原因となった呼び出し)は、それらの必要な制限の範囲外です。 – tsj

+0

これは、フィオナを使ってシェイプファイルの境界を抽出したので、「これは私にとっては意味がない」と言いました。上で、私は 'set_extent'メソッドを使って試してみました - この行を一緒に削除することを含めて - でもそれはうまくいきませんでした。 – Jhonny

+0

あなたは与えられた範囲にかかわらず、カートリストが受け入れる範囲外です。繰り返しますが、私はこれらのパッケージに精通していませんが、 'proj'の構築では、32.5Mの' false_easting'を提供しています。おそらく、これはx方向のシフトをもたらすだろう。 'set_extent'の呼び出しでは、このようなシフトによって訂正される32.5Mに近い2つの値があります。 cartopyのドキュメントを見ると、 'set_extent'への入力は(x0、x1、y0、y1)ですが、(x0、y0、x1、y1)を与えたように見えるので、' false_easting'は許容範囲に戻ってきます。 – tsj

答えて

1

私は2つの問題を発見しました。 1つはset_extentへの呼び出しの制限の間違った指定です、documentation[x0,x1,y0,y1]を指定する必要があります。[x0,y0,x1,y1]を指定したようです。 他の問題は、私が知る限り、カートリストの制限と思われます。エラーメッセージに記載されている制限外の予測は常に失敗し、その制限はハードコードされているようです。ソース(this line in their latest release)を編集するだけで、-2e7-4e7に変更することができます(上限も同様)。これらの修正の後、あなたのプロットは問題なく生成されます。 Deutschland 新しいset_extent行:

ax.set_extent([32458044.649189778*0.975, 32465287.307457082*1.025,5556418.748046352*0.9, 556415,3.5456742775*1.1],proj) 

あなたはまた、あなたのTransverseMercatorcentral_longitude=9.0を設定することもできます、あなたのシェープファイルで指定されているもののようです。

これについて開発者に連絡することをお勧めします。これらの境界を設定する理由があるか、より良い回避策があるかもしれません。また、後のリリースで境界を広げることもあります。

更新

あなたの限界もmunicipalitiesの最初に基づいて設定されているように見える:

In [34]: municipalities[0].bounds 
Out[34]: (32458044.649189778, 5556418.748046352, 32465287.307457082, 5564153.5456742775) 

しかし、他の要素が異なる境界を持ちます。すべての境界の最小値/最大値に基づいて、実際の図面に制限を適用することができます。municipalities

bnd = np.array([i.bounds for i in municipalities]) 
x0,x1 = np.min(bnd[:,0]),np.max(bnd[:,2]) 
y0,y1 = np.min(bnd[:,1]),np.max(bnd[:,3]) 
ax.set_extent([x0,x1,y0,y1],proj) 
+0

努力していただきありがとうございます。@tsj。私がフィオナを使って得たシェイプファイルの境界が、国が合う最小の正方形を私にちょうど与えないのは、私が理解していない。上記の私のコードで投影(まだ)が誤って指定されていますか?フィオナ、またはカートリに間違いはありますか?少なくとも不一致は? – Jhonny

+0

また、地球の円周は約40e6なので、カートシール側から+/- 2e7で十分であるはずです。多分、あなたの入力を変更して32.5e6の代わりに-7.5e6を参照するような方法があるかもしれません。それで、カートリストのソースを編集する必要はありません。 – tsj

+0

したがって、ジオメトリメソッドは実際にshapleyオブジェクトを返します。なぜ、私は 'bounds'を呼び出すことができますか? cartopyは '(x0、x1、y0、y1)を期待している間に'(x0、y0、x1、y1) 'を返すので、あなたがコードスニペットに表示しているときに再構成する必要があります。それを得ました、ありがとう - 魅力のように動作します。 – Jhonny

関連する問題