Eu tenho uma matriz de pontos de latitude longitude que define os limites de uma área. Eu gostaria de criar um polígono com base nesses pontos e plotar o polígono em um mapa e preenchê-lo. Atualmente, meu polígono parece consistir em muitos patches que conectam todos os pontos, mas a ordem dos pontos não está correta e quando eu tento preencher o polígono eu fico com uma área de aparência estranha (veja anexo).

Eu classifico meus pontos de longitude-latitude (array mypolyXY) de acordo com o centro do polígono, mas meu palpite é que isso não está correto:

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]))

Eu plotei os locais dos pontos (círculos pretos) e meus polígonos (manchas coloridas) usando

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)

Minha pergunta é: como posso fechar meu polígono com base em minha matriz de pontos de latitude de longitude?

ATUALIZAR: Eu testei um pouco mais sobre como plotar o polígono. Eu removi a rotina de classificação e usei os dados na ordem em que eles ocorrem no arquivo. Isso parece melhorar o resultado, mas como @tcaswell mencionou, a forma do polígono ainda se solta (veja o novo enredo). Espero que possa haver uma rotina de caminho / polígono que possa resolver meu problema e mesclar todas as formas ou caminhos dentro dos limites do polígono. Sugestões são muito bem vindas.

ATUALIZAÇÃO 2:

Agora tenho uma versão funcional do meu script que é baseada nas sugestões de @Rutger Kassies e Roland Smith. Acabei lendo o Shapefile usando org que funcionava relativamente bem. Funcionou bem para o arquivo lmes_64.shp padrão, mas quando eu usei arquivos LME mais detalhados, onde cada LME poderia consistir em vários polígonos, este script quebrou. Eu teria que encontrar uma maneira de mesclar os vários polígonos para nomes LME idênticos para fazer esse trabalho. Eu anexei o meu script que acabei no caso de alguém iria dar uma olhada nisso. Eu aprecio muito os comentários de como melhorar este script ou torná-lo mais genérico. Esse script cria os polígonos e extrai dados na região de polígono que leio em um arquivo netcdf. A grade do arquivo de entrada é -180 a 180 e -90 a 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()

Imagem final anexada. Obrigado por toda ajuda.

 Felicidades, Trond