ValueError 在 Cartopy 方位角等距地图投影中绘制 contourf

问题描述 投票：0回答：1

我想在 Cartopy 方位角等距地图投影中使用 

contourf
。对于上下文，我正在尝试绘制全球信号的传播时间（以 h 为单位）。粗略地说，这有点像我试图让我的情节看起来像：

（图片来源@htonchia，https://github.com/SciTools/cartopy/issues/1421）

但是，当我尝试绘制它时，它给出了错误： 

line 185, in geos_multipolygon_from_polygons  ValueError: Sequences of multi-polygons are not valid arguments

重现：

# Data
# Longitudes of stations
longs = [-171.7827,  -171.7827,   179.1966,   179.1966,  -159.7733,  -159.7733,
  174.7043,   174.7043,   172.9229,   172.9229,   159.9475,   159.9475,
 -157.4457, -157.4457,   146.24998,  146.24998, -169.5292,  -169.5292,
  166.652,    166.652,   -155.5326,  -155.5326,  -158.0112,  -158.0112,
 -177.3698,  -177.3698,   144.8684,   166.7572,   166.7572,   117.239,
  117.239,    125.5791,   125.5791,   110.5354,   110.5354,   144.4382,
  144.4382,   138.20406,  138.20406, -176.6842,  -176.6842,   121.4971,
  121.4971,   126.62436,  126.62436,  -64.0489,   -64.0489,  -123.3046,
 -123.3046,  -110.7847,  -110.7847,   -90.2861,   -90.2861,  -106.4572,
 -106.4572,  -106.4572,  -147.8616,  -147.8616,  -147.8616,  -104.0359,
 -104.0359,   -95.83812,  -95.83812,  -70.7005,   -70.7005,    98.9443,
   98.9443,   -88.2763,   -88.2763,   -61.9787,  -61.9787,   -78.4508,
  -78.4508,  -175.385  ]

# Latitudes of stations
lats = [-13.9085,   -13.9085,    -8.5259,    -8.5259,   -21.2125,   -21.2125,
 -41.3087,   -41.3087,     1.3549,     1.3549,    -9.4387,    -9.4387,
   2.0448,     2.0448,   -20.08765,  -20.08765,   16.7329,    16.7329,
  19.2834,    19.2834,    19.7573,    19.7573,    21.42,      21.42,
  28.2156,    28.2156,    13.5893,   -77.8492,   -77.8492,   -32.9277,
 -32.9277,     7.0697,     7.0697,   -66.2792,   -66.2792,   -89.9289,
 -89.9289,    36.54567,   36.54567,   51.8823,    51.8823,    24.9735,
  24.9735,    37.47768,   37.47768, -64.7744,   -64.7744,    44.5855,
  44.5855,    32.3098,    32.3098,    -0.6742,    -0.6742,    34.94591,
  34.94591,   34.94591,   64.873599,  64.873599,  64.873599,  44.1212,
  44.1212,    29.96478,   29.96478,  -29.011,    -29.011,     18.8141,
  18.8141,    20.2263,    20.2263,   -38.0568,   -38.0568,     0.2376,
   0.2376,   -20.57    ]

# Time (h) signal detected after eruption
travel_time_h = [ 0.95296297,  0.95332528,  1.49046297,  1.4905475,   1.67046297, 
1.67026972, 2.3705475,   2.37046297,  2.60249194,  2.60240741,  2.7537963,   2.75360306,
  3.00943639,  3.00935186,  3.65610306,  3.65601852,  3.93165861,  3.93157408,
 16.13526972,  4.43074074,  4.61268519,  4.6130475,   4.6730475,   4.67296297,
  5.01026972,  5.01046297,  5.20768519,  5.96546297,  5.9655475,   6.49693639,
  6.49685186,  6.40324074,  6.40332528,  6.53740741,  6.53721417,  7.12074074,
  7.1205475,   7.34546297,  7.34499194,  7.26157408,  7.26221417,  7.64546297,
  7.6455475,   8.13407408,  8.13388083,  7.97693639,  7.97740741,  8.05082528,
  8.05101852,  8.00240741,  8.00221417,  8.65943639,  8.65907408,  8.41907408,
  8.41776972,  8.42722222,  8.94324074,  8.9430475,   8.94333333,  9.2555475,
  9.25601852,  8.99240741,  8.99249194,  9.26851852,  9.26749194,  9.16165861,
  9.16185186,  9.41990741,  9.41999194,  9.30851852,  9.31360306,  9.82324074,
  9.82332528,  0.        ]

然后我插入数据以用于

contourf

import matplotlib.pyplot as plt
import numpy as np
from cartopy import crs as ccrs
from scipy.interpolate import griddata

# Interpolate for contour
X, Y = np.meshgrid(longs, lats)
Z = griddata((longs, lats), travel_time_h, (X, Y), method='linear')

并尝试使用

Cartopy
绘制它：

# Initialize figure
fig = plt.figure(figsize=(10, 8))
projLae = ccrs.AzimuthalEquidistant(central_longitude=-175.385, central_latitude=-20.57)
ax = plt.subplot(1, 1, 1, projection=projLae)

# Plot contour first as background
start_h, end_h, interval_h = 0.0, 10.0, 0.5
levels = np.arange(start=start_h, stop=end_h, step=interval_h)  # levels of contour
contour = ax.contourf(X, Y, Z, levels=levels, vmin=start_h, vmax=end_h, transform=ccrs.PlateCarree())

# Add colorbar for contour
cbar = fig.colorbar(contour, orientation='horizontal')   
cbar.ax.set_xlabel(f"Time [hr]")

# Plot station locations
ax.scatter(longs, lats, s=8, marker='*', color='red', transform=ccrs.PlateCarree())

# Plot map details
ax.coastlines()
ax.set_global()
plt.show()

不确定发生了什么，如果是

ax.contourf
问题和/或
Cartopy
方位角等距投影问题。我正在使用 Cartopy 版本 0.21.1.

感谢任何帮助！

python-3.x matplotlib cartopy
1个回答
0
投票

有两个主要步骤需要修改。首先，创建用于数据插值的网格网格，以及数据插值步骤。有关详细信息，请参阅代码中的注释。

更新代码：

import matplotlib.pyplot as plt
import numpy as np
from cartopy import crs as ccrs
from scipy.interpolate import griddata

# Data
# Longitudes of stations
longs = [-171.7827,  -171.7827,   179.1966,   179.1966,  -159.7733,  -159.7733,
  174.7043,   174.7043,   172.9229,   172.9229,   159.9475,   159.9475,
 -157.4457, -157.4457,   146.24998,  146.24998, -169.5292,  -169.5292,
  166.652,    166.652,   -155.5326,  -155.5326,  -158.0112,  -158.0112,
 -177.3698,  -177.3698,   144.8684,   166.7572,   166.7572,   117.239,
  117.239,    125.5791,   125.5791,   110.5354,   110.5354,   144.4382,
  144.4382,   138.20406,  138.20406, -176.6842,  -176.6842,   121.4971,
  121.4971,   126.62436,  126.62436,  -64.0489,   -64.0489,  -123.3046,
 -123.3046,  -110.7847,  -110.7847,   -90.2861,   -90.2861,  -106.4572,
 -106.4572,  -106.4572,  -147.8616,  -147.8616,  -147.8616,  -104.0359,
 -104.0359,   -95.83812,  -95.83812,  -70.7005,   -70.7005,    98.9443,
   98.9443,   -88.2763,   -88.2763,   -61.9787,  -61.9787,   -78.4508,
  -78.4508,  -175.385  ]

# Latitudes of stations
lats = [-13.9085,   -13.9085,    -8.5259,    -8.5259,   -21.2125,   -21.2125,
 -41.3087,   -41.3087,     1.3549,     1.3549,    -9.4387,    -9.4387,
   2.0448,     2.0448,   -20.08765,  -20.08765,   16.7329,    16.7329,
  19.2834,    19.2834,    19.7573,    19.7573,    21.42,      21.42,
  28.2156,    28.2156,    13.5893,   -77.8492,   -77.8492,   -32.9277,
 -32.9277,     7.0697,     7.0697,   -66.2792,   -66.2792,   -89.9289,
 -89.9289,    36.54567,   36.54567,   51.8823,    51.8823,    24.9735,
  24.9735,    37.47768,   37.47768, -64.7744,   -64.7744,    44.5855,
  44.5855,    32.3098,    32.3098,    -0.6742,    -0.6742,    34.94591,
  34.94591,   34.94591,   64.873599,  64.873599,  64.873599,  44.1212,
  44.1212,    29.96478,   29.96478,  -29.011,    -29.011,     18.8141,
  18.8141,    20.2263,    20.2263,   -38.0568,   -38.0568,     0.2376,
   0.2376,   -20.57    ]

# Time (h) signal detected after eruption
travel_time_h = [ 0.95296297,  0.95332528,  1.49046297,  1.4905475,   1.67046297, 
1.67026972, 2.3705475,   2.37046297,  2.60249194,  2.60240741,  2.7537963,   2.75360306,
  3.00943639,  3.00935186,  3.65610306,  3.65601852,  3.93165861,  3.93157408,
 16.13526972,  4.43074074,  4.61268519,  4.6130475,   4.6730475,   4.67296297,
  5.01026972,  5.01046297,  5.20768519,  5.96546297,  5.9655475,   6.49693639,
  6.49685186,  6.40324074,  6.40332528,  6.53740741,  6.53721417,  7.12074074,
  7.1205475,   7.34546297,  7.34499194,  7.26157408,  7.26221417,  7.64546297,
  7.6455475,   8.13407408,  8.13388083,  7.97693639,  7.97740741,  8.05082528,
  8.05101852,  8.00240741,  8.00221417,  8.65943639,  8.65907408,  8.41907408,
  8.41776972,  8.42722222,  8.94324074,  8.9430475,   8.94333333,  9.2555475,
  9.25601852,  8.99240741,  8.99249194,  9.26851852,  9.26749194,  9.16165861,
  9.16185186,  9.41990741,  9.41999194,  9.30851852,  9.31360306,  9.82324074,
  9.82332528,  0.        ]

# Create ranged values for making gridmesh
lon_min, lon_max, lat_min, lat_max = min(longs), max(longs), min(lats), max(lats)
grid_x, grid_y = 32, 24   # grid spacings
alons = np.linspace(lon_min, lon_max, grid_x)
alats = np.linspace(lat_min, lat_max, grid_y)

# Original bad code: X, Y = np.meshgrid(longs, lats)
X, Y = np.meshgrid(alons , alats)

# Interpolating for Z, need proper `fill_value`.
# - Available data points: longs, lats, travel_time_h
# - Places to interpolate for data: (X, Y)
# If interpolation fails, `fill_value` is used.
Z = griddata((longs, lats), travel_time_h, (X, Y), method='linear', fill_value=0.01) #linear, cubic

# Set projection
projLae = ccrs.AzimuthalEquidistant(central_longitude=-175.385, central_latitude=-20.57)

# Set figure, axes for plotting
fig = plt.figure(figsize = (10,8))
ax = fig.add_subplot(111, projection=projLae)

# Plot contour
start_h, end_h, interval_h = 0.0, 10.0, 1
levels = np.arange(start=start_h, stop=end_h, step=interval_h)  # levels of contour

contour = ax.contourf(X, Y, Z, levels=levels, vmin=start_h, vmax=end_h, transform=ccrs.PlateCarree(), zorder=20, alpha=0.75)

# Add colorbar for contour
cbar = fig.colorbar(contour, orientation='horizontal', shrink=0.45)   
cbar.ax.set_xlabel(f"Time [hr]")

# Plot station locations
ax.scatter(longs, lats, s=8, marker='*', color='red', transform=ccrs.PlateCarree(), zorder=30)

# Plot map details
ax.coastlines()
ax.set_global()

plt.show()

输出图：

不幸的是，地块在日期变更线的经度处有小的白色间隙。可以通过添加循环数据并重新生成轮廓来消除此白色间隙。相关代码和最终剧情如下

# Add cyclic data to eliminate the white gap
import cartopy.util as cutil
cdata, clon, clat = cutil.add_cyclic(Z, X, Y)
contour = ax.contourf(clon, clat, cdata, levels=levels, vmin=start_h, vmax=end_h, transform=ccrs.PlateCarree(), zorder=20, alpha=0.75)

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.