Mam tablicę punktów długości i szerokości geograficznej, która określa granice obszaru. Chciałbym utworzyć wielokąt na podstawie tych punktów i narysować wielokąt na mapie i wypełnić go. Obecnie mój wielokąt wydaje się składać z wielu łat, które łączą wszystkie punkty, ale kolejność punktów nie jest prawidłowa i kiedy próbuję wypełnić wielokąt, otrzymuję dziwnie wyglądający obszar (patrz załączony).

Sortuję punkty długości i szerokości geograficznej (tablica mypolyXY) zgodnie ze środkiem wielokąta, ale przypuszczam, że nie jest to poprawne:

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

Korzystam z wykreślania lokalizacji punktów (czarne kółka) i moich wielokątów (kolorowe łaty)

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)

Moje pytanie brzmi: jak mogę zamknąć swój wielokąt na podstawie mojej tablicy punktów długości i szerokości geograficznej?

AKTUALIZACJA: Sprawdziłem trochę więcej o tym, jak narysować wielokąt. Usunąłem procedurę sortowania i po prostu użyłem danych w kolejności, w jakiej występują w pliku. Wydaje się, że poprawia to wynik, ale jak już wspominałem @tcaswell, kształt wielokąta nadal się podcina (patrz nowy wykres). Mam nadzieję, że może istnieć procedura ścieżki / wielokąta, która może rozwiązać mój problem i połączyć wszystkie kształty lub ścieżki w granicach wielokąta. Sugestie są bardzo mile widziane.

AKTUALIZACJA 2:

Mam teraz działającą wersję mojego skryptu opartą na sugestiach @Rutger Kassies i Roland Smith. Skończyło się na tym, że czytałem Shapefile używając org, który działał stosunkowo dobrze. Sprawdziło się to w przypadku standardowego pliku lmes_64.shp, ale kiedy użyłem bardziej szczegółowych plików LME, gdzie każdy LME mógł składać się z kilku wielokątów, ten skrypt się zepsuł. Musiałbym znaleźć sposób na połączenie różnych wielokątów dla identycznych nazw LME, aby to działało. Dołączam swój skrypt, na który skończyłem, gdyby ktoś mógł go obejrzeć. Bardzo cenię komentarze dotyczące ulepszenia tego skryptu lub uczynienia go bardziej ogólnym. Ten skrypt tworzy wielokąty i wyodrębnia dane w obszarze wielokąta, który czytam z pliku netcdf. Siatka pliku wejściowego wynosi od -180 do 180 i od -90 do 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()

Dołączony obraz końcowy. Dziękuję za wszelką pomoc.

 Pozdrawiam, Trond