在python中将时间序列数据转换为网格化3D数组

Question

我有一个包含多年时间序列数据的数据框，并且每天包含不同纬度位置的变量值。对于给定的一天，变量记录在不同的位置。以下是我在 python pandas 中读取的数据帧的 sn-p：

               lat      lon         variable  
Date                                                            
2017-12-31  12.93025  59.9239     10.459373     
2019-12-31  12.53044  43.9229     12.730064     
2019-02-28  12.37841  33.9245     37.487683  

我想：

1. 将其网格化为 2x2.5 度分辨率
2. 制作一个包含网格数据及其时间变化的 3D 数组。 我想得到一个网格数据集作为具有形状（时间，纬度，经度）的数组。这是因为我以特定分辨率网格化的数据帧必须与全球气象数据进行比较分辨率为 2x2.5 度。 （另外，我的数据集不会记录所有日期的所有位置的数据，并且必须在创建最终数组时处理丢失的数据）。

我研究了 geopandas、xarray 和 histogram2d 来对数据进行网格化。我还使用 histigram2d 函数成功地对数据进行了网格化。但是，只能实现缺少时间信息的 2D 阵列，这使我的分析成为一项挑战。我知道，理想情况下，我应该将时间维度连接到我的 2D 数组，但考虑到并非所有位置都始终记录数据，因此很难做到这一点。

这就是我使用 histogram2d 函数创建 1 度网格单元的方式：

**

#Plot histogram2d - for gridding the data:
df=df_in['2019'] #taking one year at a time
# Test data, globally distributed
lat_r = df['lat']
lon_r = df['lon']
z_r = df['variable']
lat = np.array(lat_r)
lon = np.array(lon_r)
z = np.array(z_r)
    
# Create binning
binlon = np.linspace(-180,180, 361)
binlat = np.linspace(-90, 90, 181)
zz, xx, yy = np.histogram2d(lon, lat, bins=(binlon, binlat), weights=z, normed=False)
counts, _, _= np.histogram2d(lon, lat, bins=(binlon, binlat))\

# Workaround for zero count values tto not get an error.
# Where counts == 0, zi = 0, else zi = zz/counts
zi = np.zeros_like(zz)
zi[counts.astype(bool)] = zz[counts.astype(bool)]/counts[counts.astype(bool)]
zi = np.ma.masked_equal(zi, 0)

#Final, gridded data:
hist = zi.T # shape(180,360)

**

在这方面的任何帮助将不胜感激。

Answer 1

我最终制作了示例数据并处理了 2D 和 3D 案例。我将从您已经使用的 2D 案例开始，因为 3D 案例的扩展非常简单。

二维

首先，让我们创建一些随机样本数据。请注意，我在这里导入了以后需要的所有内容

import numpy as np
import matplotlib.pyplot as plt
import cartopy
from cartopy.crs import PlateCarree
from matplotlib.colors import Normalize

def create2Ddata():
    '''Makes some random data'''

    N = 2000
    lat = 10 * np.random.rand(N) + 40
    lon = 25 * np.random.rand(N) - 80
    z = np.sin(4*np.pi*lat/180.0*np.pi) + np.cos(8*np.pi*lon/180.0*np.pi)

    return lat, lon, z

# Create Data
lat, lon, z = create2Ddata()

这将用作我们想要使用直方图函数绘制的一些随机、分散的地理空间数据。下一步是创建有意义的分箱，然后进行实际分箱。

def make2dhist(lon, lat, z, latbins, lonbins):
    '''Takes the inputs and creates 2D histogram'''
    zz, _, _ = np.histogram2d(lon, lat, bins=(
        lonbins, latbins), weights=z, normed=False)
    counts, _, _ = np.histogram2d(lon, lat, bins=(lonbins, latbins))\

    # Workaround for zero count values to not divide by zero.
    # Where counts == 0, zi = 0, else zi = zz/counts
    zi = np.zeros_like(zz)
    zi[counts.astype(bool)] = zz[counts.astype(bool)] / \
        counts[counts.astype(bool)]
    zi = np.ma.masked_equal(zi, 0)

    return lonbins, latbins, zi

# Make bins
latbins = np.linspace(np.min(lat), np.max(lat), 75)
lonbins = np.linspace(np.min(lon), np.max(lon), 75)

# Bin the data
_, _, zi = make2dhist(lon, lat, z, latbins, lonbins)

然后，我们将分散数据和分箱数据绘制如下。

def plotmap():
    '''background map plotting'''

    ax = plt.gca()
    ax.add_feature(cartopy.feature.LAND, zorder=0, edgecolor='None',
                   linewidth=0.5, facecolor=(0.8, 0.8, 0.8))
    ax.spines['geo'].set_linewidth(0.75)


fig = plt.figure()

# Just plot the scattered data
ax = plt.subplot(211, projection=PlateCarree())
plotmap()
plt.scatter(lon, lat, s=7, c=z, cmap='rainbow')

# Plot the binned 2D data
ax = plt.subplot(212, projection=PlateCarree())
plotmap()
plt.pcolormesh(
    lonbins, latbins, zi.T, shading='auto', transform=PlateCarree(),
    cmap='rainbow')
plt.show()

Figure 2D, not allowed to embed figures yet...

顶部是分散的数据，底部是分箱数据。

---

3D

让我们继续 3D 案例。同样，让我们​​创建一些随时间变化的随机分散数据：

def create3Ddata():
    ''' Make random 3D data '''
    N = 8000
    lat = 10 * np.random.rand(N) + 40
    lon = 25 * np.random.rand(N) - 80
    t = 10 * np.random.rand(N)

    # Linearly changes sign of the cos+sin wavefield
    z = (t/5 - 1) * (np.sin(2*2*np.pi*lat/180.0*np.pi)
                     + np.cos(4*2*np.pi*lon/180.0*np.pi))

    return lat, lon, t, z


# Create Data
lat, lon, t, z = create3Ddata()

现在，这里不再使用histogram2d，而是使用histogramdd，这只是同一函数的N维版本。

def make3dhist(lon, lat, t, z, latbins, lonbins, tbins):
    '''Takes the inputs and creates 3D histogram just as the 2D histogram
    function'''
    zz, _ = np.histogramdd(
        np.vstack((lon, lat, t)).T,
        bins=(lonbins, latbins, tbins),
        weights=z, normed=False)

    counts, _ = np.histogramdd(
        np.vstack((lon, lat, t)).T,
        bins=(lonbins, latbins, tbins))
    # Workaround for zero count values tto not get an error.
    # Where counts == 0, zi = 0, else zi = zz/counts
    zi = np.zeros_like(zz)
    zi[counts.astype(bool)] = zz[counts.astype(bool)] / \
        counts[counts.astype(bool)]
    zi = np.ma.masked_equal(zi, 0)
    return lonbins, latbins, tbins, zi

# Create bins
latbins = np.linspace(np.min(lat), np.max(lat), 75)
lonbins = np.linspace(np.min(lon), np.max(lon), 75)
tbins = np.linspace(np.min(t), np.max(t), 5)

# Bin the data
_, _, _, zi = make3dhist(lon, lat, t, z, latbins, lonbins, tbins)

最后，我们在各自的时间箱中并排绘制分散数据和分箱数据。请注意用于确保容易观察到时间变化的标准化。 请注意，共有三个循环（我可以将它们放在一个循环中，但这样更便于阅读）。

1. 第一个循环及时对数据进行分箱，并将分箱后的数据分别绘制在一个切片中。
2. 第二个循环使用之前的 2D 直方图函数将数据按时间分箱，然后在空间分箱，并为每个时间箱绘制一个切片。
3. 第三个函数使用上面已经 3D 分箱的数据，并通过访问 3D 矩阵中的切片来绘制切片。

# Normalize the colors so that variations in time are easily seen
norm = Normalize(vmin=-1.0, vmax=1.0)

fig = plt.figure(figsize=(12, 10))

# The scattered data in time bins
# Left column
for i in range(4):
    ax = plt.subplot(4, 3, 3*i + 1, projection=PlateCarree())
    plotmap()

    # Find points in time bins
    pos = np.where((tbins[i] < t) & (t < tbins[i+1]))

    # Plot scatter points
    plt.title(f'{tbins[i]:0.2f} < t < {tbins[i+1]:0.2f}')
    plt.scatter(lon[pos], lat[pos], c=z[pos], s=7, cmap='rainbow', norm=norm)
    plt.colorbar(orientation='horizontal', pad=0.0)

# Center column
for i in range(4):
    ax = plt.subplot(4, 3, 3*i + 2, projection=PlateCarree())
    plotmap()
    plt.title(f'{tbins[i]:0.2f} < t < {tbins[i+1]:0.2f}')

    # Find data points in time bins
    pos = np.where((tbins[i] < t) & (t <= tbins[i+1]))

    # Bin the data for each time bin separately
    _, _, zt = make2dhist(lon[pos], lat[pos], z[pos], latbins, lonbins)
    plt.pcolormesh(
        lonbins, latbins, zt.T, shading='auto', transform=PlateCarree(),
        cmap='rainbow', norm=norm)
    plt.colorbar(orientation='horizontal', pad=0.0)

# Right column
for i in range(4):
    ax = plt.subplot(4, 3, 3*i + 3, projection=PlateCarree())
    plotmap()
    plt.title(f'{tbins[i]:0.2f} < t < {tbins[i+1]:0.2f}')
    plt.pcolormesh(
        lonbins, latbins, zi[:, :, i].T, shading='auto', transform=PlateCarree(),
        cmap='rainbow', norm=norm)
    plt.colorbar(orientation='horizontal', pad=0.0)

plt.show()

Figure 3D, not allowed to embed figures yet...

在左列中，是分散的、随机的地理空间数据，其中标题表示分箱。在中心列中，使用“手动”时间分档数据的 2D 直方图。在右栏中，使用 3D 直方图分箱的切片。 正如预期的那样，中心列和右列显示完全相同的内容。

希望这能解决您的问题。