如何用matplotlib画dem图像

https://github.com/rveciana/geoexamples
import matplotlib.pyplot as plt

from mpl_toolkits.basemap import Basemap

def hillshade(array, azimuth, angle_altitude):  # 刻画视觉3D效果

    x, y = gradient(array)

    slope = pi/2. - arctan(sqrt(x*x + y*y))

    aspect = arctan2(-x, y)

    azimuthrad = azimuth*pi / 180.

    altituderad = angle_altitude*pi / 180.

    shaded = sin(altituderad) * sin(slope)\

     + cos(altituderad) * cos(slope)\

     * cos(azimuthrad - aspect)

    return 255*(shaded + 1)/2

fig,ax=plt.subplots()

m=Basemap(projection='merc',llcrnrlat=29,llcrnrlon=102,urcrnrlat=30,urcrnrlon=104,ax=ax)

#ds = gdal.Open('/home/loong1/data/Asia_topo30.grd')

ds = gdal.Open('leshan.grd') #事先准备好grd数据文件

band = ds.GetRasterBand(1)

arr = band.ReadAsArray()#可以从grd文件里读取部分内容，但这里读取了全部

lon,lat=103,29.5

x,y=m(lon,lat)

hs_array = hillshade(arr,315, 45)

m.scatter(x,y,marker='o',c='r')

(x1,y1)=m(102,29)

(x2,y2)=m(104,30)

m.imshow(hs_array,cmap='Greys',origin='upper',extent=[x1,x2,y1,y2]) # extent参数指定图片的坐标范围，叠加到了basemap给定的投影底图上

m.drawparallels(arange(29,30,0.5),labels=[0,1,1,0])

m.drawmeridians(arange(102,104,1),labels=[1,0,0,1])

plt.show()