排序多邊形座標繪圖

Question

我有一個模型網格組成的許多單元格，我想繪製一個陰影多邊形的matplotlibbasemap。

使用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) 

當試圖這樣做時，我得到了一個十字交叉的模式，我認爲它必須對我提供給多邊形函數的座標進行排序。

我認爲這與我如何提取外部座標有關，或者只是在創建最終的多邊形時繪製座標列表的順序。

如果出現問題，是否有巧妙的排序方法？

繪製的多邊形

特寫

Answer 1

SOOOO ...顯然shapely.geometry.Polygon方法繪製與所有內部網格多邊形座標，其中我意識到由於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)) 

Answer 2

我同意存在的網格座標一些misarrangement。 grid_lons是如何創建的？可能使用Pyproj和Shapely幾何體的更簡潔的方法是使用相對較新的功能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) 

座標的完整性應該通過轉換來保存（假設它們是開始的）。