从边界点创建封闭多边形

Question

我有一个经度 - 纬度点阵列，它定义了一个区域的边界。我想创建一个基于这些点的多边形，并在地图上绘制多边形并填充它。目前，我的多边形似乎由许多连接所有点的补丁组成，但点的顺序不正确，当我尝试填充多边形时，会看到一个奇怪的外观区域（请参阅附件）。

我排序根据多边形的中心我的经纬度点（mypolyXY阵列），但我的猜测是，这是不正确的：

cent=(np.sum([p[0] for p in mypolyXY])/len(mypolyXY),np.sum([p[1] for p in mypolyXY])/len(mypolyXY)) 
# sort by polar angle 
mypolyXY.sort(key=lambda p: math.atan2(p[1]-cent[1],p[0]-cent[0])) 

我绘制的点位置（黑圈）和我的多边形（色斑）使用

scatter([p[0] for p in mypolyXY],[p[1] for p in mypolyXY],2) 
p = Polygon(mypolyXY,facecolor=colors,edgecolor='none') 
ax.add_artist(p) 

我的问题是：我怎么能收我的多边形根据我的经纬度点阵列上？

更新： 我测试了一些关于如何绘制多边形。我删除了排序例程，并按照它们在文件中出现的顺序使用数据。这似乎改善了结果，但正如@tcaswell所提到的那样，多边形的形状本身仍然是凹陷的（参见新图）。我希望可以有一个路径/多边形例程，可以解决我的问题，并在多边形的边界内合并所有形状或路径。建议是非常受欢迎的。

更新2：

我现在我的脚本是基于@Rutger Kassies和罗兰·史密斯建议的工作版本。我最终使用org读取了Shapefile，这个工作比较好。它对标准的lmes_64.shp文件运行良好，但是当我使用更详细的LME文件，其中每个LME可以由多个多边形组成时，该脚本就会崩溃。我将不得不找到一种方法来合并不同的多边形以获得相同的LME名称来完成这项工作。我附上了我的剧本，最后我会介绍一下，以防有人看到它。我非常感谢有关如何改进此脚本或使其更通用的评论。该脚本创建多边形并提取我从netcdf文件读取的多边形区域内的数据。输入文件的网格为-180至180和-90至90.

import numpy as np 
import math 
from pylab import * 
import matplotlib.patches as patches 
import string, os, sys 
import datetime, types 
from netCDF4 import Dataset 
import matplotlib.nxutils as nx 
from mpl_toolkits.basemap import Basemap 
import ogr 
import matplotlib.path as mpath 
import matplotlib.patches as patches 


def getLMEpolygon(coordinatefile,mymap,index,first): 

    ds = ogr.Open(coordinatefile) 
    lyr = ds.GetLayer(0) 
    numberOfPolygons=lyr.GetFeatureCount() 

    if first is False: 
     ft = lyr.GetFeature(index) 
     print "Found polygon:", ft.items()['LME_NAME'] 
     geom = ft.GetGeometryRef() 

     codes = [] 
     all_x = [] 
     all_y = [] 
     all_XY= [] 

     if (geom.GetGeometryType() == ogr.wkbPolygon): 
      for i in range(geom.GetGeometryCount()): 

      r = geom.GetGeometryRef(i) 
      x = [r.GetX(j) for j in range(r.GetPointCount())] 
      y = [r.GetY(j) for j in range(r.GetPointCount())] 

      codes += [mpath.Path.MOVETO] + (len(x)-1)*[mpath.Path.LINETO] 
      all_x += x 
      all_y += y 
      all_XY +=mymap(x,y) 


     if len(all_XY)==0: 
      all_XY=None 
      mypoly=None 
     else: 
      mypoly=np.empty((len(all_XY[:][0]),2)) 
      mypoly[:,0]=all_XY[:][0] 
      mypoly[:,1]=all_XY[:][3] 
    else: 
     print "Will extract data for %s polygons"%(numberOfPolygons) 
     mypoly=None 
    first=False 
    return mypoly, first, numberOfPolygons 


def openCMIP5file(CMIP5name,myvar,mymap): 
    if os.path.exists(CMIP5name): 
     myfile=Dataset(CMIP5name) 
     print "Opened CMIP5 file: %s"%(CMIP5name) 
    else: 
     print "Could not find CMIP5 input file %s : abort"%(CMIP5name) 
     sys.exit() 
    mydata=np.squeeze(myfile.variables[myvar][-1,:,:]) - 273.15 
    lonCMIP5=np.squeeze(myfile.variables["lon"][:]) 
    latCMIP5=np.squeeze(myfile.variables["lat"][:]) 

    lons,lats=np.meshgrid(lonCMIP5,latCMIP5) 

    lons=lons.flatten() 
    lats=lats.flatten() 
    mygrid=np.empty((len(lats),2)) 
    mymapgrid=np.empty((len(lats),2)) 

    for i in xrange(len(lats)): 
     mygrid[i,0]=lons[i] 
     mygrid[i,1]=lats[i] 
     X,Y=mymap(lons[i],lats[i]) 
     mymapgrid[i,0]=X 
     mymapgrid[i,1]=Y 

    return mydata, mygrid, mymapgrid 

def drawMap(NUM_COLORS): 

    ax = plt.subplot(111) 
    cm = plt.get_cmap('RdBu') 
    ax.set_color_cycle([cm(1.*j/NUM_COLORS) for j in range(NUM_COLORS)]) 
    mymap = Basemap(resolution='l',projection='robin',lon_0=0) 

    mymap.drawcountries() 

    mymap.drawcoastlines() 
    mymap.fillcontinents(color='grey',lake_color='white') 
    mymap.drawparallels(np.arange(-90.,120.,30.)) 
    mymap.drawmeridians(np.arange(0.,360.,60.)) 
    mymap.drawmapboundary(fill_color='white') 

    return ax, mymap, cm 

"""Edit the correct names below:""" 

LMEcoordinatefile='ShapefileBoundaries/lmes_64.shp' 
CMIP5file='tos_Omon_CCSM4_rcp85_r1i1p1_200601-210012_regrid.nc' 

mydebug=False 
doPoints=False 
first=True 


"""initialize the map:""" 
mymap=None 
mypolyXY, first, numberOfPolygons = getLMEpolygon(LMEcoordinatefile, mymap, 0,first) 
NUM_COLORS=numberOfPolygons 
ax, mymap, cm = drawMap(NUM_COLORS) 


"""Get the CMIP5 data together with the grid""" 
SST,mygrid, mymapgrid = openCMIP5file(CMIP5file,"tos",mymap) 

"""For each LME of interest create a polygon of coordinates defining the boundaries""" 
for counter in xrange(numberOfPolygons-1): 

    mypolyXY,first,numberOfPolygons = getLMEpolygon(LMEcoordinatefile, mymap,counter,first) 

    if mypolyXY != None: 
     """Find the indices inside the grid that are within the polygon""" 
     insideBoolean = plt.mlab.inside_poly(np.c_[mymapgrid[:,0],mymapgrid[:,1]],np.c_[mypolyXY[:,0],mypolyXY[:,1]]) 
     SST=SST.flatten() 
     SST=np.ma.masked_where(SST>50,SST) 

     mymapgrid=np.c_[mymapgrid[:,0],mymapgrid[:,1]] 
     myaverageSST=np.mean(SST[insideBoolean]) 

     mycolor=cm(myaverageSST/SST.max()) 
     scaled_z = (myaverageSST - SST.min())/SST.ptp() 
     colors = plt.cm.coolwarm(scaled_z) 

     scatter([p[0] for p in mypolyXY],[p[1] for p in mypolyXY],2) 

     p = Polygon(mypolyXY,facecolor=colors,edgecolor='none') 
     ax.add_artist(p) 

     if doPoints is True: 

      for point in xrange(len(insideBoolean)): 
       pointX=mymapgrid[insideBoolean[point],0] 
       pointY=mymapgrid[insideBoolean[point],1] 
       ax.scatter(pointX,pointY,8,color=colors) 
       ax.hold(True) 


if doPoints is True: 
    colorbar() 
print "Extracted average values for %s LMEs"%(numberOfPolygons) 
plt.savefig('LMEs.png',dpi=300) 
plt.show() 

附加的最终图像。感谢所有帮助。

干杯，特朗德

Answer 1

我推荐使用原始Shapefile，它的格式适合存储多边形。至于OGR替代你可以使用身材匀称，或多边形输出到高铁总站等

import ogr 
import matplotlib.path as mpath 
import matplotlib.patches as patches 
import matplotlib.pyplot as plt 

ds = ogr.Open('lmes_64.shp') 
lyr = ds.GetLayer(0) 
ft = lyr.GetFeature(38) 
geom = ft.GetGeometryRef() 
ds = None 

codes = [] 
all_x = [] 
all_y = [] 

if (geom.GetGeometryType() == ogr.wkbPolygon): 
    for i in range(geom.GetGeometryCount()): 

    r = geom.GetGeometryRef(i) 
    x = [r.GetX(j) for j in range(r.GetPointCount())] 
    y = [r.GetY(j) for j in range(r.GetPointCount())] 

    codes += [mpath.Path.MOVETO] + (len(x)-1)*[mpath.Path.LINETO] 
    all_x += x 
    all_y += y 

if (geom.GetGeometryType() == ogr.wkbMultiPolygon): 
    codes = [] 
    for i in range(geom.GetGeometryCount()): 
    # Read ring geometry and create path 
    r = geom.GetGeometryRef(i) 
    for part in r: 
     x = [part.GetX(j) for j in range(part.GetPointCount())] 
     y = [part.GetY(j) for j in range(part.GetPointCount())] 
     # skip boundary between individual rings 
     codes += [mpath.Path.MOVETO] + (len(x)-1)*[mpath.Path.LINETO] 
     all_x += x 
     all_y += y 

carib_path = mpath.Path(np.column_stack((all_x,all_y)), codes)  
carib_patch = patches.PathPatch(carib_path, facecolor='orange', lw=2) 

poly1 = patches.Polygon([[-80,20],[-75,20],[-75,15],[-80,15],[-80,20]], zorder=5, fc='none', lw=3) 
poly2 = patches.Polygon([[-65,25],[-60,25],[-60,20],[-65,20],[-65,25]], zorder=5, fc='none', lw=3) 


fig, ax = plt.subplots(1,1) 

for poly in [poly1, poly2]: 
    if carib_path.intersects_path(poly.get_path()): 
     poly.set_edgecolor('g') 
    else: 
     poly.set_edgecolor('r') 

    ax.add_patch(poly) 

ax.add_patch(carib_patch) 
ax.autoscale_view() 

也会检出Fiona（包装为OGR），如果你想很容易的Shapefile处理。

Answer 2

具有点的数组是不够的。您需要知道点的订单。通常，多边形的点顺序为。所以你从第一点到第二点画一条线，然后从第二点到第三点画一条线。

如果您的列表不是按顺序排列，则需要额外的信息才能制作顺序列表。

shapefile（请参阅documentation）包含形状列表，如空形状，点，多边形，多边形，其变体还包含Z和M（测量）坐标。所以只是倾销积分不会。你必须把它们分成不同的形状，并渲染你感兴趣的。在这种情况下可能是PolyLine或Polygon。有关这些形状的数据格式，请参阅上面的链接。请记住，文件的某些部分是big-endian，而另一些则是little-endian。真是一团糟。

我会建议使用struct模块来解析二进制.shp文件，因为再根据该文件，一个多边形的点为了和他们形成闭合链（最后一点是相同的作为第一）。

您可以尝试使用当前坐标列表的另一件事是从一个点开始，然后在列表中寻找相同的点。这些相同点之间的所有内容应该是一个多边形。这可能不是万无一失的，但看看它有多远。

Answer 3

有一个博客张贴在这里它谈论shape文件和底图：http://www.geophysique.be/2013/02/12/matplotlib-basemap-tutorial-10-shapefiles-unleached-continued/

如果你是热衷于尝试一下，cartopy也可能是一种选择。从shapefile中绘制数据被设计为相当简单：

import matplotlib.pyplot as plt 

import cartopy.crs as ccrs 
import cartopy.io.shapereader as shpreader 

# pick a shapefile - Cartopy makes it easy to access Natural Earth 
# shapefiles, but it could be anything 
shapename = 'admin_1_states_provinces_lakes_shp' 
states_shp = shpreader.natural_earth(resolution='110m', 
             category='cultural', name=shapename) 

。

# states_shp is just a filename to a shapefile 
>>> print states_shp 
/Users/pelson/.local/share/cartopy/shapefiles/natural_earth/cultural/110m_admin_1_states_provinces_lakes_shp.shp 

。

# Create the mpl axes of a PlateCarree map 
ax = plt.axes(projection=ccrs.PlateCarree()) 

# Read the shapes from the shapefile into a list of shapely geometries. 
geoms = shpreader.Reader(states_shp).geometries() 

# Add the shapes in the shapefile to the axes 
ax.add_geometries(geoms, ccrs.PlateCarree(), 
        facecolor='coral', edgecolor='black') 

plt.show()