如何以非均匀分辨率在物理距离坐标中绘制数组

Question

问题

我想处理模拟 64*30（XZ 平面）多云大气的一些数据。

具体来说，我想显示一些物理属性的图像，例如消光，作为位置的函数（物理距离单位（km），而不是“x点数，z点数”/“像素”坐标，如果不清楚，将进一步举例说明).

• 我有水平分辨率：55m
• 我知道垂直轴从 0 到 1000m。垂直分辨率不是恒定的，它在云层内加倍（因为我们想要更高的精度）。
• 即：z ∈ [0,400]m，分辨率为50m，即从索引z[0]到z[9] - 在云下
  z ∈ [400,850]m，分辨率为 25m，即从索引 z[9] 到 z[26] - 在云中
  z ∈ [850,1000]m，分辨率又是 50m，即从索引 z[26] 到 z[29] - 在云上

问题基本上是在不规则网格中显示数组。我不知道最有效的方法是什么。
有什么方法是我没有想到的吗？

我做了什么
我的数据组织为：

X1 Z1 EXTINCTION ... other parameters ...
X1 Z2 EXTINCTION ... other parameters ...
 .
 .
 .
X1 Z30 EXTINCTION ...other parameters ...
X2 Z1 EXTINCTION ...other parameters ...
 .
 .
 .
X64 Z30 EXTINCTION ...other parameters ...

我已经将每一行的消光值存储在某个数组

EXT

中，以及每行的x和z坐标

X

＆

Z

接下来，我将

EXT

的值添加到表示空间坐标的 64*30 矩阵中。 我在范围循环中使用 2 来遵循数据中行索引的顺序，它由 k 表示。

Extinction = np.zeros(64,30)
k = 0 #count the line
for x in range (0,64): 
        for z in range(0,30) : 
            Extinction [x,z] = EXT[k]  #Add the value of extinction for the corresponding x,z position     
            k += 1 #next line

最后，我用

matplotlib.pyplot.imshow

显示结果

fig, ax = plt.subplots()
im = ax.imshow(Extinction.T,cmap='gnuplot',origin='lower') 

因为我希望我的轴以米为单位显示距离，知道分辨率和像素数，我只使用：

extent(55*64,0,1000)

作为

ax.imshow()

中的论点

当然，结果图仅适用于恒定的垂直分辨率。 Extent() 只均匀地拉伸图。由于分辨率从 50m 变为 25m，高度将不正确。。

“蛮力”法

我尝试采用看起来 easy 的方式。我不确定它是否正确。
我通过人为地将云外的点数加倍（实质上是将分辨率加倍）来修改我的

Extinction

数组。我这样做是为了让垂直分辨率在任何地方都保持不变。我的想法是，既然消光系数在空气中总是0，那应该不是问题。我做了以下事情：

import numpy as np
import matplotlib.pyplot as plt

Alebdo = np.zeros(64,30)
k = 0
for x in range (0,64): 
        for z in range(0,17): #below the cloud
            Extinction[x,z] = 0 
        k += 9 #we doubled each line. The line number corresponding the beginning of the cloud is 9     
        
        for z in range(17,36): #cloudy zone
            Extinction[x,z] = EXT[k]     
            k += 1
        
        for z in range(36,40): #above the cloud
            Extinction[x,z] = 0     
        k += 2

这似乎可行，但是：

1. 云层有点“挤”
2. 不知道是不是物理上的“正确”
3. 一些属性（不像反照率或消光）在云下不是恒定的。所以我不能只把云层从空气中分离出来，然后在其他地方设置一个常数值。例如，云层下方的辐射场将不均匀。
4. 这不是很实用，因为我必须根据我使用的模拟手动更改很多参数。

插值

我之前没怎么用interpolate。但据我所知，它在这种不规则网格的情况下很强大：

from scipy.interpolate import griddata

xg = np.linspace(0,64,41)
zg = np.arange(0,30,1)
xgrid, zgrid = np.meshgrid(xg, zg)
EXT_interpolate = griddata((X, Z), EXT, (xgrid, zgrid), method='nearest')

im=ax.imshow(xgrid, zgrid, EXT_interpolate, cmap='gnuplot')
f=fig.colorbar(im, ax=ax,shrink=0.8)

我的问题：

1. 我的图表以“数组索引”坐标显示数据。如何显示以米为单位的距离？我不能像以前那样简单地使用

   extent()

   ，

   contourf

   有什么等价物吗？
2. 我应该在哪里插入我的数据？使用 41 个垂直点是有意义的（因为这是我在分辨率为 25m 时拥有的点数：40*25m = 1000m 并添加水平 z=0）。
   
   编辑： 感谢 Paul Brodersen 的评论，我选择了插值：

#X was previously a list of index : X=[1,2,3....,64] 
#I changed it to distances in km since I know the resolution is 55m X=[0,0.055,0.110 ,...]
#Same for Z, which is now Z=[0,0.050,0.0100,0.0150, ..., 0.350,0.400,0.425,0.450,...]
# So it now accounts for the change in resolution

from scipy.interpolate import griddata

#Grid on which we will interpolate
xg = np.linspace(0,3.465,0.055) #don't touch the x resolution
zg = np.arange(0,1,100)#change z resolution so that it's uniform, 10m here.
xgrid, zgrid = np.meshgrid(xg, zg)
EXT_interpolate = griddata((X, Z), EXT, (xgrid, zgrid), method='nearest')

#Plot the interpolation
im=ax.imshow(xgrid, zgrid, EXT_interpolate, cmap='gnuplot2')
f=fig.colorbar(im, ax=ax,shrink=0.8)

#I still have to add extent(0,3.465,0,1) if I want to display distances and not the index of my grid points

Answer 1

引用pcolormesh文档

创建一个带有不规则矩形网格的伪彩色图。

下面的代码比您需要的要长，因为我想向您展示着色选项如何改变轴限制的解释方式。

另外，对于

shading='flat'

，你必须传递一个“更小”的颜色矩阵，我只是删除了最后一行，最后一列，你会执行某种平均/插值。

最终，网格线显示网格是不规则的，在生产中你可能会删除

ec='w', lw=0.1, antialiased=1

参数

from matplotlib.pyplot import show, subplots
from numpy import arange, concatenate, cos, meshgrid

x = 55.0*arange(64)
z = concatenate((
    0.00 + 50*arange( 8),
    400. + 25*arange(18),
    850. + 50*arange( 4),
))
X, Z = meshgrid(x, z)
E = cos((X/3000-Z/1200)) # unimportant

fig, ((none, flat), (nearest, gouraud)) = subplots(2, 2, figsize=(11,4.5), layout='constrained')

none.pcolormesh(X, Z, E, shading=None, ec='w', lw=0.2, antialiased=1)
none.set_aspect(1)
none.set_title('No shading')

flat.pcolormesh(X, Z, E[:-1,:-1], shading='flat', ec='w', lw=0.2, antialiased=1)
flat.set_aspect(1)
flat.set_title('Flat shading')

nearest.pcolormesh(X, Z, E, shading='nearest', ec='w', lw=0.2, antialiased=1)
nearest.set_aspect(1)
nearest.set_title('Nearest shading')

gouraud.pcolormesh(X, Z, E, shading='gouraud', ec='w', lw=0.2, antialiased=1)
gouraud.set_aspect(1)
gouraud.set_title('Gouraud shading')

show()