排序多边形坐标绘图

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) 

坐标的完整性应该通过转换来保存（假设它们是开始的）。