2

多くのセルで構成されるモデル グリッドがあり、シェーディングされたポリゴンをmatplotlib basemap.

を使用してpyproj、最初にポイントを投影してから、shapely.geometryのPolygonクラスを使用してポリゴンを作成し、グリッドの外部座標を抽出しました。次に、プロット関数に渡すためにそれらを WGS84 に戻します。

grid_x_mesh, grid_y_mesh = pyproj.transform(wgs84, nplaea, grid_lons, grid_lats)
grid_x = grid_x_mesh.ravel()
grid_y = grid_y_mesh.ravel()
grid_poly = Polygon(zip(grid_x, grid_y))
grid_x, grid_y = grid_poly.exterior.coords.xy
grid_plons, grid_plats = pyproj.transform(nplaea, wgs84, grid_x, grid_y)

次に、matplotlib.basemapメソッドを使用して、WSG84 座標を地図投影 (この場合は nplaea) に投影し、

grid_poly_x, grid_poly_y = m(grid_plons, grid_plats)
grid_poly_xy = zip(grid_poly_x, grid_poly_y)
grid_poly = Polygon(grid_poly_xy, facecolor='red', alpha=0.4)
plt.gca().add_patch(grid_poly)

そうしようとすると、十字形のパターンが得られます。これは、ポリゴン関数に指定した座標の順序付けを行う必要があると想定しています。

これは、外部座標を抽出した方法、またはプロットする最終的なポリゴンを作成したときの座標リストの順序に関係していると思います。

それが問題である場合、これらを適切に注文する賢い方法はありますか?

プロットされたポリゴン ここに画像の説明を入力

閉じる ここに画像の説明を入力

4

2 に答える 2

1

グリッド座標の配置に誤りがあることに同意します。どのようにgrid_lons作成されましたか？おそらく、Shapely ジオメトリで Pyproj を使用するよりクリーンな方法は、比較的新しい関数を使用することshapely.ops.transformです。例えば：

import pyproj
from shapely.geometry import Polygon, Point
from shapely.ops import transform
from functools import partial

project = partial(
    pyproj.transform,
    pyproj.Proj(init='epsg:4326'),  # WGS84 geographic
    pyproj.Proj(init='epsg:3575'))  # North Pole LAEA Europe

# Example grid cell, in long/lat
poly_g = Polygon(((5, 52), (5, 60), (15, 60), (15, 52), (5, 52)))

# Transform to projected system
poly_p = transform(project, poly_g)

座標の健全性は、変換を通じて保持する必要があります (最初から健全であると仮定して)。

于 2014-02-11T22:05:45.280 に答える
0

すっごく...どうやらshapely.geometry.Polygonメソッドは、すべての内部グリッド座標でポリゴンを描画していたようです。これは、 edメッシュ座標配列grid_plonsとgrid_plats同じ長さであるために実現しました。np.ravel()

メソッドに渡す前に、メッシュ座標配列から外部座標を手動で抽出するだけで済みましたPolygon（以下を参照）。ただし、これを行うには、よりきれいで一般的な方法があると思います。

手動抽出方法:

grid_x_mesh, grid_y_mesh = pyproj.transform(wgs84, nplaea, grid_lons, grid_lats)

# The coordinates must be ordered in the order they are to be drawn
[grid_x.append(i) for i in grid_x_mesh[0,:]]
[grid_x.append(i) for i in grid_x_mesh[1:-1,-1]]
# Note that these two sides of the polygon are appended in reverse
[grid_x.append(i) for i in (grid_x_mesh[-1,:])[::-1]]
[grid_x.append(i) for i in (grid_x_mesh[1:-1,0])[::-1]]

[grid_y.append(i) for i in grid_y_mesh[0,:]]
[grid_y.append(i) for i in grid_y_mesh[1:-1,-1]]
[grid_y.append(i) for i in (grid_y_mesh[-1,:])[::-1]]
[grid_y.append(i) for i in (grid_y_mesh[1:-1,0])[::-1]]

grid_poly = Polygon(zip(grid_x, grid_y))
于 2014-03-13T14:56:47.280 に答える