如何在 python 中使用 cartopy 跨越而不跨越日期变更线绘制多边形？

Question

我有几个多边形要一起可视化。一些多边形穿过日期变更线。我不知道如何将多边形一起可视化。

我已经尝试了用于创建形状多边形的点的多种排列方式，尝试了 central_longitude 投影参数（例如 cartopy set extent with central_longitude=180）和经度的转换（它们本身以度为单位东，0：360）。后者似乎影响最大。也就是说，当我不进行转换时，太平洋多边形是正确的，但海湾多边形没有显示出来（图Correct Pacific, no Gulf。关闭校正，两个多边形都出现了，但太平洋多边形不正确（图 Incorrect Pacific）。感谢您提供的所有帮助。

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

Answer 1

要获得预期效果，必须在代码的所有部分正确使用 CRS。通常，我们的数据编码在一个 CRS 中，通常是地理数据（经度、纬度）。下面代码中，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()

Answer 2

在我之前的回答中，多边形不是根据大地测量法生成并绘制在地图上的。所以，需要这个答案。请注意，需要一些技巧才能获得良好的结果。

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

剧情：

如何在 python 中使用 cartopy 跨越而不跨越日期变更线绘制多边形？

How to plot polygons across, and not across, dateline with cartopy in python?

python

gis

matplotlib

shapely

cartopy