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

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, форма многоугольника все еще подрывается (см. Новый график). Я надеюсь, что может быть процедура path/polygon, которая могла бы решить мою проблему и объединить все формы или пути в границах многоугольника. Предложения приветствуются.

UPDATE 2:

У меня сейчас рабочая версия моего сценария, который основан на предложениях @Rutger Kassies и Роланд Смит. Я закончил читать 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() 

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

Cheers, Тронд

Answer 1

Я рекомендую использовать оригинальный Shapefile, который находится в формате, подходящем для хранения полигонов. В качестве альтернативы OGR можно использовать стройные, или экспортировать многоугольник WKT и т.д.

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

Имея массив точек не достаточно. Вы должны знать заказ пунктов. Обычно точки многоугольника задаются последовательно. Таким образом, вы рисуете линию от первой точки ко второй, затем от второй до третьей и т. Д.

Если ваш список не в последовательном порядке, вам нужна дополнительная информация, чтобы иметь возможность составлять последовательный список.

Файл формы (см. documentation) содержит список фигур, таких как нулевая фигура, точка, полилиния, многоугольник, варианты, содержащие также координаты Z и M (измерение). Так что просто сбросить очки не удастся. Вы должны разделить их в разных формах и отображать те, которые вас интересуют. В этом случае, вероятно, PolyLine или Polygon. См. Ссылку выше для формата данных для этих фигур.Имейте в виду, что некоторые части файла являются большими, а другие - маленькими. Какой беспорядок.

Я предложил бы использовать struct модуль для разбора двоичного .shp файла, потому что снова в соответствии с документацией, точки одного многоугольника являются в порядке, и они образуют замкнутую цепь (последний пункт такой же как первый).

Еще одна вещь, которую вы могли бы попробовать с вашим текущим списком координат - начать с точки, а затем искать ту же точку дальше в списке. Все между этими точками должно быть одним полигоном. Это, вероятно, не является надежным, но посмотрите, как далеко он вам достанется.

Answer 3

Существует сообщение в блоге здесь, который говорит о шейпфайлах и базовой карте: http://www.geophysique.be/2013/02/12/matplotlib-basemap-tutorial-10-shapefiles-unleached-continued/

Если вы заинтересованы, чтобы попробовать его, cartopy также может быть вариантом. Построение данных из шейп-файла было достаточно простым:

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