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