python - 如何使用 ESRI ASCII 栅格格式文件(如 Matlab 中的 mapshow)在 Python 中绘制栅格 map ？

Question

例如，我有一个 ESRI ASCII 栅格格式的 DEM 文件，如下所示:

ncols 480
nrows 450
xllcorner 378923
yllcorner 4072345
cellsize 30
nodata_value -32768
43 2 45 7 3 56 2 5 23 65 34 6 32 54 57 34 2 2 54 6 
35 45 65 34 2 6 78 4 2 6 89 3 2 7 45 23 5 8 4 1 62 ...

我想绘制一个栅格 map 来显示地形。我知道它可以通过 Matlab 中的 mapshow 来实现，比如

[Z,R] = arcgridread(filename);
mapshow(Z,R,'DisplayType','Surface')

但是如何在 Python 中实现呢？以及如果坐标系是英国国家网格，是否可以在Python中在栅格 map 上添加一个shapefile图层(例如多边形文件)？

Answer 1

嗯，虽然有两个人认为这个问题很有用，但没人愿意回答这个问题。所以，我现在可以自己回答了。 首先，读取 ESRI ASCII 栅格文件 ('file.asc')

import numpy as np
# read header
file_name = 'file.asc'
header_rows = 6 # six rows for header information
header = {} # store header information including ncols, nrows, xllcorner, yllcorner, cellsize, NODATA_value
row_ite = 1
with open(file_name, 'rt') as file_h:
     for line in file_h:
        if row_ite <= header_rows:
             line = line.split(" ", 1)
             header[line[0]] = float(line[1])
        else:
             break
        row_ite = row_ite+1
# read data array
data_array = np.loadtxt(file_name, skiprows=header_rows, dtype='float64')

然后计算出DEM网格的范围，即它的四个边的坐标。

left = header['xllcorner']
right = header['xllcorner']+header['ncols']*header['cellsize']
bottom = header['yllcorner']
top = header['yllcorner']+header['nrows']*header['cellsize']
map_extent = (left, right, bottom, top)

然后你可以用 DEM 数组和它的范围绘制 map 。

import matplotlib.pyplot as plt
fig, ax = plt.subplots(1)
img = plt.imshow(data_array, extent=map_extent)

如果要添加颜色条，可能需要以下代码:

from mpl_toolkits.axes_grid1 import make_axes_locatable
divider = make_axes_locatable(ax)
cax = divider.append_axes(loc='right', size='3%', pad=0.05,)
cbar = plt.colorbar(img, cax=cax)

使用这些代码，您还可以在 Python 中定义自己的 arcgridread 和 mapshow 函数。