Python：用 np.array 数据覆盖 shapefile

Question

我有一个美国的 shapefile，我有一个 m x n 笛卡尔数据数组，代表每个像素的温度。我能够加载 shapefile 并绘制它：

import shapefile as shp
import matplotlib.pyplot as plt

sf = shp.Reader("/path/to/USA.shp")

plt.figure()
for shape in sf.shapeRecords():
    for i in range(len(shape.shape.parts)):
        i_start = shape.shape.parts[i]
        if i==len(shape.shape.parts)-1:
            i_end = len(shape.shape.points)
        else:
            i_end = shape.shape.parts[i+1]
        x = [i[0] for i in shape.shape.points[i_start:i_end]]
        y = [i[1] for i in shape.shape.points[i_start:i_end]]
        plt.plot(x,y, color = 'black')
plt.show()

而且我能够读取我的数据并绘制它：

import pickle
from matplotlib import pyplot as mp
Tfile = '/path/to/file.pkl'
with open(Tfile) as f:
    reshapeT = pickle.load(f)
mp.matshow(reshapeT)

问题是 reshapeT 的尺寸为 536 x 592，并且是美国的子域。但是，我确实有关于 reshapeT 网格左上角的信息（纬度/经度）以及每个像素之间的间距 (0.01)

我的问题是：如何将 reshapeT 数据叠加到 shapefile 域之上？

Answer 1

如果我对你的理解是正确的，你想在绘制的 shapefile 的特定部分上覆盖一个 536x592 numpy 数组。我建议您使用带有 extent 参数的 Matplotlib 的 imwshow() 方法，它允许您将图像放置在绘图中。

你绘制 shapefile 的方式很好，但是，如果你有可能使用 geopandas，它将大大简化事情。绘制 shapefile 将减少为以下几行：

import geopandas as gpd
sf = gpd.read_file("/path/to/USA.shp")
ax1 = sf.plot(edgecolor='black', facecolor='none')

正如您之前所做的那样，现在让我们加载数组数据：

import pickle
Tfile = '/path/to/file.pkl'
with open(Tfile) as f:
    reshapeT = pickle.load(f)

现在为了能够在正确的位置绘制 numpy 数组，我们首先需要计算它的范围（它将覆盖的区域以坐标表示）。你提到你有关于 top-left 角和分辨率 (0.01) 的信息——这就是我们所需要的。在下文中，我假设有关 top-left 角的 lat/lon 信息保存在 top_left_lat 和 top_left_lon 变量中。范围需要在一个元组中传递，每个边都有一个值（按左、右、下、上的顺序）。

因此，我们的范围可以计算如下：

extent_mat = (top_left_lon, top_left_lon + reshapeT.shape[1] * 0.01, top_left_lat - reshapeT.shape[0] * 0.01, top_left_lat)

最后，我们将矩阵绘制到同一个坐标轴对象上，ax1，我们已经将形状文件绘制到计算范围：

# Let's turn off autoscale first. This prevents
# the view of the plot to be limited to the image
# dimensions (instead of the entire shapefile). If you prefer
# that behaviour, just remove the following line
ax1.autoscale(False)

# Finally, let's plot!
ax1.imshow(reshapeT, extent=extent_mat)

Python：用 np.array 数据覆盖 shapefile

Python: Overlaying shapefile with np.array data

python

matplotlib

shapefile