[FIXED] How to draw raster map in Python with a ESRI ASCII Raster format file like mapshow in Matlab?

Issue

For example, I have a DEM file in ESRI ASCII Raster format like this:

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 ...

I want to draw a raster map to display the topography. I know it can be achieved via mapshow in Matlab like

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

But how to do it in Python? And if the coordinate system is British National Grid, is it possible to add a shapefile layer (such as a polygon file) on the raster map in Python?

Solution

Well, nobody would like to answer this question although two people think the question is useful. So, I can answer it by myself now. Firstly, read the ESRI ASCII Raster file ('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')

Then calculate the extent of the DEM grid, namely the coordinates of its four edges.

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)

Then you can draw a map with DEM array and its extent.

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

If you want to add a colorbar, you may need the following codes:

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)

With these codes, you can also define your own arcgridread and mapshow functions in Python.

Answered By - Sheldon