Создать замкнутый многоугольник из граничных точек

У меня есть массив точек долготы-широты, который определяет границы области. Я хотел бы создать многоугольник на основе этих точек, нанести его на карту и заполнить его. В настоящее время мой многоугольник, кажется, состоит из множества патчей, которые соединяют все точки, но порядок точек неправильный, и когда я пытаюсь заполнить многоугольник, я получаю странную область (см. Приложение).

Я сортирую свои точки долготы и широты (массив 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 и Roland Smith. Я закончил читать Shapefile, используя org, который работал относительно хорошо. Это работало хорошо для стандартного файла 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()

Окончательное изображение прилагается. Спасибо за помощь.

 Ура, Тронд

Ответы на вопрос(3)

Ваш ответ на вопрос