题目：

data.txt 中的数据表示地球表面各个经纬度坐标的海拔高度，其中每行三个浮点数，分别是经度[-180, 180]，纬度[-90, 90]和高度。数据中经纬度的间隔是 0.1 度，高度的单位是米。请设计合适的颜色映射函数，实现地形图的经纬度投影。

使用python matplotlib库中的basemap toolkit

Anaconda源和pip源都安装不成功，需要在这里下载，本地安装

另外，matplotlib不像matlab中有专门适应地形图的色盘demcmap，因此需要自己设置颜色映射，即0以下数值映射到浅蓝到深蓝代表海平面以下，0以上采用terrain色盘，采用经纬度投影（'cyl’）, contourf()函数绘制

head(data.txt)

世界地形图（经纬度投影）

在第 1 题的基础上，实现地形图的方位角等距投影，圆锥等面积（亚尔勃斯）投影和圆柱等角度（墨卡托）投影。

设置投影参数为方位角等距投影(‘aeqd’)，圆锥等面积（亚尔勃斯）投影(‘aea’)和圆柱等角度（墨卡托）投影(‘merc’)即可

圆锥等面积（亚尔勃斯）投影

方位角等距投影

圆柱等角度（墨卡托）投影

分别在上一题的三种投影方式所生成的可视化结果中添加交互：使用鼠标选择地球上两个位置 A、B，绘制两地之间的飞行航线，假设飞行航线始终为测地线（即球面最短距离，由大圆所确定的球面距离，你需要计算这条路径）。注意：由于飞机在平流层飞行，你不需要考虑地表海拔对飞机的阻隔。

定义一个array存放每次点击事件下的地图坐标，每两次点击，将地图坐标转化为经纬度坐标，使用drawgreatcircle()绘制测地线，并清空array

鼠标点击绘制测地线

测试环境：Python-3.6，basemap-1.1.0

code

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors

class FixPointNormalize(matplotlib.colors.Normalize):

    def __init__(self, vmin=None, vmax=None, sealevel=0, col_val = 0.21875, clip=False):
        # sealevel is the fix point of the colormap (in data units)
        self.sealevel = sealevel
        # col_val is the color value in the range [0,1] that should represent the sealevel.
        self.col_val = col_val
        matplotlib.colors.Normalize.__init__(self, vmin, vmax, clip)

    def __call__(self, value, clip=None):
        x, y = [self.vmin, self.sealevel, self.vmax], [0, self.col_val, 1]
        return np.ma.masked_array(np.interp(value, x, y))

colors_undersea = plt.cm.terrain(np.linspace(0, 0.17, 56))
colors_land = plt.cm.terrain(np.linspace(0.25, 1, 200))
colors = np.vstack((colors_undersea, colors_land))
cut_terrain_map = matplotlib.colors.LinearSegmentedColormap.from_list('cut_terrain', colors)

data = np.loadtxt("data.txt",delimiter=' ')
mlon = np.copy(data[:,0]).reshape(1801,3601)
mlat = np.copy(data[:,1]).reshape(1801,3601)
mdem = np.copy(data[:,2]).reshape(1801,3601)

norm = FixPointNormalize(sealevel=0, vmax=8500)

m = Basemap(projection='merc',llcrnrlat=-80,urcrnrlat=80, llcrnrlon=-180,urcrnrlon=180,lat_ts=20,resolution='c')
lat, lon = m(mlon,mlat)
m.drawcoastlines()
m.contourf(lat, lon, mdem,cmap = cut_terrain_map, norm = norm, levels = list(range(-11000, 8500, 250)))
m.drawparallels(np.arange(-90,90,30),labels=[1,0,0,0])
m.drawmeridians(np.arange(-180,180,30),labels=[0,0,0,1])
cbar = m.colorbar()
plt.title('world elevation map')

linepoints = np.array([])

def onclick(event):
    if event.button == 3:
        global linepoints
        x = event.xdata
        y = event.ydata
        xx, yy = m(x, y, inverse=True)
        linepoints = np.append(linepoints, xx)
        linepoints = np.append(linepoints, yy)
        if np.size(linepoints) == 4:
            m.drawgreatcircle(lon1 = linepoints[0], lat1 = linepoints[1], lon2=linepoints[2], lat2=linepoints[3], color = 'red',linewidth = 3)
            linepoints = np.array([])
            plt.show()

plt.gcf().canvas.mpl_connect("button_press_event", onclick)
plt.show()