2

我有一个由许多单元格组成的模型网格，我想在matplotlib basemap.

使用，我首先投影点，然后使用'类pyproj创建多边形以从中提取网格的外部坐标。然后我将它们恢复为 WGS84 以传递给我的绘图函数：shapely.geometryPolygon

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方法是使用所有内部网格坐标绘制多边形，我意识到这是由于grid_plons和具有与'ed 网格坐标数组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 回答