python-3.x Cartopy方位等距Map投影等值线绘制的数值误差

nhjlsmyf  于 2023-02-26  发布在  Python
关注(0)|答案(1)|浏览(264)

我想在Cartopy Azimuth EquidistanceMap投影中使用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。
我很感激你的帮助!

xlpyo6sf

xlpyo6sf1#

有两个主要步骤需要修改。首先,创建用于数据插值的gridmesh,以及数据插值步骤。详细信息请参见代码中的注解。
更新代码:

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 columns, rows
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)

相关问题