Как построить многоугольники поперек, а не поперек линии дат с картографией в python? - PullRequest
Новая
       143

Как построить многоугольники поперек, а не поперек линии дат с картографией в python?

0 голосов
/ 04 августа 2020

У меня есть несколько полигонов, которые я пытаюсь визуализировать вместе. Некоторые многоугольники пересекают линию дат. Я не могу понять, как визуализировать многоугольники вместе.

Я пробовал несколько вариантов расположения точек, используемых для создания фигурного многоугольника, экспериментировал с параметром проекции central_longitude (например, экстент набора карт с central_longitude = 180 ) и преобразования долгот (которые изначально находятся в градусах восточной долготы, 0: 360). Последнее, кажется, имеет наибольший эффект. То есть, когда я не выполняю преобразование, многоугольник Pacifi c правильный, но многоугольник Gulf не отображается (рис. Correct Pacifi c, no Gulf . Отключение исправление, оба полигона отображаются, но один из них Pacifi c неверен (рис. Incorrect Pacifi c). Спасибо за любую помощь.

MWE 

import cartopy.crs as ccrs
from shapely.geometry import Point, Polygon
plt.style.use('seaborn-dark-palette')

## degrees lon (0 : 360), lat (-90 : 90)
polys = {'SPNA': [(250, 25), (280, 25), (302, 10)],
           'EP': [(178, 48), (227, 48), (227, 24)]}
           
crs       = ccrs.PlateCarree(central_longitude=180)
fig, axes = plt.subplots(figsize=(14,8), subplot_kw={'projection': crs})
axes.stock_img()

for loc, poly in polys.items():
    pts = []; lons, lats = [], []
    for lon, lat in poly:
        ## uncomment this to produce the difference in the two attached images
        # lon = lon - 360 if lon >= 180 else lon # 0:360 to -180:180
        pt  = Point(lon, lat)
        pts.append(pt)
        lons.append(pt.x)
        lats.append(pt.y)
    shp     = Polygon(pts)
    print (shp)

    axes.scatter(lons, lats, transform=ccrs.PlateCarree())
    axes.add_geometries([shp], crs=ccrs.PlateCarree())
plt.show()

Ответы [ 2 ]

2 голосов
/ 04 августа 2020

Чтобы получить то, что вы ожидаете, необходимо правильно использовать CRS во всех частях вашего кода. Обычно наши данные кодируются в одной CRS, обычно это geographi c (долгота, широта в градусах). В приведенном ниже коде crs0 - это CRS данных, а crs180 - это CRS проекции карты. Поскольку два CRS различны, преобразование координат должно выполняться на определенных этапах, чтобы получить правильный результат. 

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from shapely.geometry import Point, Polygon

plt.style.use('seaborn-dark-palette')

# define CRS's 
crs0 = ccrs.PlateCarree(central_longitude=0)       #for coding data
lon_0 = 180  #OK for 180, for 0, you may be surprised, but it's still OK
crs180 = ccrs.PlateCarree(central_longitude=lon_0) #for plotting map

# (coding data)
# data is in `crs0`
polys = {'SPNA': [(250, 25), (280, 25), (302, 10)],
           'EP': [(178, 48), (227, 48), (227, 24)]}

# (specs of plotting, use `crs180`)
fig, axes = plt.subplots(figsize=(8,5), subplot_kw={'projection': crs180})
axes.stock_img()

for loc, poly in polys.items():
    pts = []
    lons, lats = [], []
    for lon, lat in poly:
        pt  = Point(lon, lat)
        # transform data (x,y) to what we need for plotting
        px, py = crs180.transform_point(lon, lat, crs0)
        pts.append( [px, py] )
        lons.append(px)
        lats.append(py)

    shp = Polygon(pts)
    #print (shp)

    # Note the options: `transform`, and `crs`
    # Our data are already in `crs180`, same as axes' CRS
    #  `transform` can be ignored in some places
    axes.scatter(lons, lats, transform=crs180)
    axes.add_geometries([shp], crs=crs180, fc="gold", ec="red")

plt.show()

обновленная карта

0 голосов
/ 04 августа 2020

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

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from shapely.geometry import Point, Polygon

# -360, is used to shift longitude of some polygons to get it work
# Note: sequence of vertices goes clock-wise for all polygons
lonshift = -360  # multiple not OK, positive not OK
polys = {'SPNA': [(250+lonshift, 25), (280+lonshift, 25), (302+lonshift, 10)],
           'EP': [(158+lonshift, 48), (227+lonshift, 48), (227+lonshift, 24)],
           'EXTRA': [(60, -40),(100, 25),(150, 10)]}  # No shift is allowed for this polygon west of CM

crs       = ccrs.PlateCarree(central_longitude=180)
fig, axes = plt.subplots(figsize=(8,5), subplot_kw={'projection': crs})
axes.stock_img()

for loc, poly in polys.items():
    pts = []
    lons, lats = [], []
    # must reverse sequence of vertices here
    poly.reverse()   # reverse to get counter clockwise sequences (in place)
    for lon, lat in poly:
        pt  = Point(lon, lat)
        pts.append(pt)
        lons.append(pt.x)
        lats.append(pt.y)

    shp = Polygon(pts)
    #print (shp)

    axes.scatter(lons, lats, transform=ccrs.PlateCarree())

    # note the use of .as_geodetic() to get geodesic as perimeter of polygons
    axes.add_geometries([shp], crs=ccrs.PlateCarree().as_geodetic(), fc="green", ec="red", alpha=0.65)
    # Filled color will be on the outside of the polygon if their vertices go clockwise in this step

plt.show()

Сюжет:

введите описание изображения здесь

...