netCDF 多边形是固定大小的光栅图像

我正在尝试从 netCDF 文件转换数据；数据包含多边形和甲烷浓度值。

from matplotlib.patches import Polygon
import cartopy.crs as ccrs
import cartopy.feature as cfeature
# %% Plotting the polygons
fig, ax = plt.subplots(figsize=(10, 10), subplot_kw={'projection': ccrs.PlateCarree()})
ax.set_extent([-9, -5, 32, 34], crs=ccrs.PlateCarree())

# Add a base map
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle=':')

# Normalize the ch4 values for the colormap
norm = mcolors.Normalize(vmin=1850, vmax=1880)
cmap = matplotlib.colormaps['rainbow']

# Plot the grid cells and color them based on ch4 values
for lat_corners, lon_corners, ch4 in zip(lat_corners, lon_corners, ch4_values):
    color = cmap(norm(ch4))
    poly = Polygon(list(zip(lon_corners, lat_corners)), facecolor=color,
                   edgecolor='grey', linewidth=0.5, alpha=0.7)
    ax.add_patch(poly)

# Add a colorbar
sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
sm.set_array([])
cbar = plt.colorbar(sm, ax=ax, orientation='horizontal', fraction=0.046, pad=0.04)
cbar.set_label('Methane Concentration (ppb)')
plt.show()

该图像非常具体，我想应用它来生成形状为 [N,M,32x32] 像素图像的数据集，其中 N 是图像数量M 包含不同形式的数据（甲烷值、风等）。

我遇到的一个大问题是将多边形之间的空间（不存在的地方）视为一个像素。

I对任何信息缺失表示歉意；这对我来说是非常新的事情。

def project_to_image(lat_corners, lon_corners, ch4_values):
    # Create a 32x32 empty image
    #image = np.zeros((32, 32, 3), dtype=np.uint8)
    #image = np.full((32, 32, 3), 255, dtype=np.uint8)
    norm = mcolors.Normalize(vmin=1850, vmax=1880)
    cmap = matplotlib.colormaps['rainbow']
    #min_value_color = cmap(norm(1850))  # Color for the minimum value in the range
    #min_value_color = (np.array(min_value_color[:3]) * 255).astype(np.uint8)  # Convert to RGB
    min_value_color = 255
    # Create a 32x32 empty image with the background color set to the minimum value color
    image = np.full((32, 32, 3), min_value_color, dtype=np.uint8)

    # Create a Basemap instance
    m = Basemap(projection='cyl', llcrnrlat=31.5, urcrnrlat=34.0, llcrnrlon=-9.0, urcrnrlon=-5.0, resolution='h')

    for lat_corner, lon_corner, ch4_value in zip(lat_corners, lon_corners, ch4_values):
        if np.any(~lon_corner.mask) and np.any(~lat_corner.mask):
            color = cmap(norm(ch4_value))
            x, y = m(lon_corner[~lon_corner.mask], lat_corner[~lat_corner.mask])

            # Limit coordinates to Basemap range
            x = np.clip(x, m.llcrnrx, m.urcrnrx)
            y = np.clip(y, m.llcrnry, m.urcrnry)


            # Adjust coordinates to fit 32x32 image
            # Map coordinates to 0-31 range
            x = (31*(x - m.llcrnrx) / (m.urcrnrx - m.llcrnrx) ).astype(int)
            y = (31*(y - m.llcrnry) / (m.urcrnry - m.llcrnry) ).astype(int)
            # Filter out of range coordinates
            #valid_idx = (x >= 0) & (x < 32) & (y >= 0) & (y < 32)
            #x = x[valid_idx]
            #y = y[valid_idx]

            # Draw polygon
            for i in range(len(x)):
                image[int(y[i]), int(x[i])] = (np.array(color[:3]) * 255).astype(np.uint8)

    return image

我正在尝试做的事情：

我试图从中复制数据的研究文章可以在这里找到：

研究文章

好的，我了解想实现的目标。想将 netCDF 文件中表示为多边形的数据转换为固定大小的 32x32 像素图像，类似于在研究论文中提供的示例。面临的挑战是正确地将多边形光栅化为图像，包括填充多边形之间的空间，并表示背景值。

以下是一种使用 Python 和适当的库（如 rasterio 用于地理空间栅格处理和 matplotlib 用于颜色图）实现此目的的方法：

import netCDF4 as nc
import numpy as np
import rasterio
from rasterio import features
from rasterio.transform import from_bounds
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize
from matplotlib import cm

# 1. 读取 netCDF 数据
data_path = "your_netcdf_file.nc"  # 替换为的 netCDF 文件路径
dataset = nc.Dataset(data_path)

# 假设的变量名为“ch4”，并且有“lat_corners”和“lon_corners”变量
ch4_values = dataset.variables['ch4'][:]
lat_corners = dataset.variables['lat_corners'][:]
lon_corners = dataset.variables['lon_corners'][:]

# 2. 设置图像参数
image_width = 32
image_height = 32
# 获取数据的范围以创建转换
min_lon = np.min(lon_corners)
max_lon = np.max(lon_corners)
min_lat = np.min(lat_corners)
max_lat = np.max(lat_corners)
transform = from_bounds(min_lon, min_lat, max_lon, max_lat, image_width, image_height)

# 3. 创建用于光栅化的函数
def rasterize_polygons(shapes, values, transform, width, height):
    """将多边形光栅化为具有给定变换和维度的栅格。

    Args:
        shapes: 几何图形列表。
        values: 与形状对应的值列表。
        transform: 用于光栅化的仿射变换。
        width: 输出栅格的宽度。
        height: 输出栅格的高度。

    Returns:
        光栅化多边形的 NumPy 数组。
    """
    raster = features.rasterize(
        ((geom, value) for geom, value in zip(shapes, values)),
        out_shape=(height, width),
        transform=transform,
        fill_value=np.nan,  # 将多边形外部的值设置为 NaN
        dtype=np.float32
    )
    return raster

# 4. 循环浏览每个时间步长（或拥有的任何其他维度）并创建图像
images = []
for i in range(len(ch4_values)):  # 假设的数据具有时间维度
    # 为当前时间步长创建多边形
    shapes = []
    for j in range(len(ch4_values[i])):
        polygon = Polygon(zip(lon_corners[i, j], lat_corners[i, j]))
        shapes.append(polygon)

    # 将多边形光栅化为图像
    raster = rasterize_polygons(shapes, ch4_values[i], transform, image_width, image_height)

    # 将 NaN 值替换为背景值或插值
    # 可以使用 np.nanmean() 或任何其他插值技术
    raster[np.isnan(raster)] = np.nanmean(ch4_values)  

    # 5. 对图像应用颜色图（可选）
    norm = Normalize(vmin=1850, vmax=1880)
    cmap = cm.get_cmap('rainbow')
    # 将颜色图应用于数据
    colored_raster = cmap(norm(raster))

    # 将颜色图数组转换为 0-255 范围内的 uint8
    image = (colored_raster[:, :, :3] * 255).astype(np.uint8)

    images.append(image)

# 现在有一个图像列表（每个时间步长一个）
# 可以将它们堆叠到一个 NumPy 数组中，或者根据需要单独处理它们
images = np.array(images)

在此代码中：

1. 我们首先从 netCDF 文件中加载数据。
2. 我们定义输出图像的参数，包括宽度、高度和仿射变换。仿射变换对于将经纬度坐标映射到图像像素很重要。
3. rasterize_polygons 函数采用几何图形、值、变换、宽度和高度作为输入，并返回一个 NumPy 数组，该数组表示光栅化的多边形。
4. 我们循环浏览每个时间步长，创建多边形，将它们光栅化为图像，然后将 NaN 值替换为选定的背景值。
5. 我们使用选定的颜色图和归一化将光栅化数据转换为彩色图像。

此代码片段应该可以帮助生成所需的数据集，其中多边形被准确地光栅化为固定大小的图像，包括填充多边形之间的空间。请记住根据的特定文件、变量名和所需背景处理调整代码。