如何在等高线图上绘制具有条件的数组？

小能豆

如何在等高线图上绘制具有条件的数组？

我使用以下代码绘制了 GPP 的全球地图：

（’lon’ 和 ‘lat’ 都是 netCDF4 属性，形状分别为 (144, ) 和 (90, )，而 ‘gpp_avg’ 是一个 numpy 数组，形状为 (90, 144) ）

import numpy as np
import netCDF4 as n4
import matplotlib.pyplot as plt

import cartopy as cart
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from mpl_toolkits.basemap import Basemap

>> gpp_avg = n4.Dataset('decadal_gpp.nc', 'r')
>> lon = gpp_avg.variables['lon'] # 144 grid cells every 2.5 degrees (east-west)
>> lat = gpp_avg.variables['lat'] # 90 grid cells every 2 degrees (north-south)

>> # Plotting data on a map with Cartopy
>> plt.figure()
>> ax = plt.axes(projection=ccrs.PlateCarree())
>> ax.coastlines() # Adding coastlines
>> ax.add_feature(cart.feature.OCEAN, zorder=100, edgecolor='k')
>> cs = ax.contourf(lon[:], lat[:], gpp_avg[:], cmap = 'Spectral')

>> cbar = plt.colorbar(cs, ax=ax) # Additional necessary information
>> cbar.set_label('g[C]/m^2/day')
>> gridl = ax.gridlines(color="black", linestyle="dotted", 
           draw_labels=True) # Adding axis labels - latitude & longitude
>> gridl.xformatter=LONGITUDE_FORMATTER
>> gridl.yformatter=LATITUDE_FORMATTER
>> gridl.xlabels_top = False
>> gridl.ylabels_right = False
>> plt.show()

我有一个 numpy 数组“ci_95_gpp”，其形状为 (90, 144)，其中包含全局地图每个网格单元的 p 值。我想在全球轮廓图上绘制 p 值大于 2 的点。

我该怎么做呢？非常感谢。

阅读 17

2024-12-11

共1个答案

小能豆

要在你已经绘制的全球地图上添加满足条件（p值 > 2）的点，可以通过在现有地图上使用 scatter 函数绘制这些点。scatter 函数可以将给定的经纬度位置上的点标记出来。以下是如何在你的代码中实现这个功能。

假设 ci_95_gpp 数组包含每个网格单元的 p 值，并且你想在这些 p 值大于 2 的位置上绘制点。

修改代码来绘制 `p值 > 2` 的点

你可以首先找到 ci_95_gpp > 2 的经纬度坐标，然后使用 ax.scatter() 绘制这些点。

修改后的代码：

import numpy as np
import netCDF4 as n4
import matplotlib.pyplot as plt
import cartopy as cart
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

# 读取数据
gpp_avg = n4.Dataset('decadal_gpp.nc', 'r')
lon = gpp_avg.variables['lon'][:]  # 经度 (144)
lat = gpp_avg.variables['lat'][:]  # 纬度 (90)
gpp_avg_data = gpp_avg.variables['gpp_avg'][:]  # GPP数据 (90, 144)

# 假设 ci_95_gpp 是已知的 90x144 数组，其中存储每个网格单元的 p 值
ci_95_gpp = np.random.rand(90, 144) * 5  # 生成一个示例 p 值数组 (你应该使用实际数据)

# 找到 p值 > 2 的网格位置
mask = ci_95_gpp > 2

# 获取经度和纬度的索引
lon_points = lon[np.where(mask)]
lat_points = lat[np.where(mask)]

# 绘制地图
plt.figure()
ax = plt.axes(projection=ccrs.PlateCarree())
ax.coastlines()  # 绘制海岸线
ax.add_feature(cart.feature.OCEAN, zorder=100, edgecolor='k')

# 绘制 GPP 数据
cs = ax.contourf(lon, lat, gpp_avg_data, cmap='Spectral')

# 添加颜色条
cbar = plt.colorbar(cs, ax=ax)
cbar.set_label('g[C]/m^2/day')

# 添加网格线
gridl = ax.gridlines(color="black", linestyle="dotted", draw_labels=True)
gridl.xformatter = LONGITUDE_FORMATTER
gridl.yformatter = LATITUDE_FORMATTER
gridl.xlabels_top = False
gridl.ylabels_right = False

# 绘制 p值 > 2 的点
ax.scatter(lon_points, lat_points, color='red', marker='o', s=10, label="p > 2")

# 显示图例
ax.legend()

# 展示结果
plt.show()

解释

数据读取：从 NetCDF 文件中读取 lon、lat 和 gpp_avg 数据。ci_95_gpp 是一个包含 p 值的数组，形状为 (90, 144)，你需要根据实际数据来替换 ci_95_gpp。
条件筛选：使用 mask = ci_95_gpp > 2 来选择那些 p 值大于 2 的网格单元。
获取坐标：np.where(mask) 返回满足条件的经纬度索引，lon_points 和 lat_points 提取出这些经纬度坐标。
绘制地图：使用 Cartopy 绘制 GPP 数据，并在这些满足 p > 2 条件的位置上使用 scatter 绘制红色圆点。
图例：通过 ax.legend() 来添加图例，标注红点是表示 p > 2 的点。

结果

运行以上代码后，你将看到在 GPP 全球地图上，p > 2 的位置会以红色小圆点的形式标注出来。如果你有多个阈值或需要自定义颜色、大小等参数，可以根据需要调整 scatter() 函数的参数。

希望这个方法对你有所帮助！

2024-12-11