.png%3Ftable%3Dblock%26id%3D99172ae2-ee44-4671-9354-b020520c82ba%26cache%3Dv2&w=3840&q=75)
.png%3Ftable%3Dblock%26id%3Da50dae53-2423-4b78-8a48-65f1e7d52926%26cache%3Dv2&w=3840&q=75)
绘制研究区域DEM地形图
1. 下载数据
全局1分钟网格中以ASCII XYZ格式提取DEM地形数据并绘图。下载数据保存为csv
2. 数据预览
3. 代码
import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import griddata import matplotlib.colors import pandas as pd from matplotlib import rcParams config = {"font.family":'Times New Roman', "font.size": 16, "mathtext.fontset":'stix'} rcParams.update(config) class FixPointNormalize(matplotlib.colors.Normalize): """ Inspired by https://stackoverflow.com/questions/20144529/shifted-colorbar-matplotlib Subclassing Normalize to obtain a colormap with a fixpoint somewhere in the middle of the colormap. This may be useful for a `terrain` map, to set the "sea level" to a color in the blue/turquise range. """ 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)) # Combine the lower and upper range of the terrain colormap with a gap in the middle # to let the coastline appear more prominently. # inspired by https://stackoverflow.com/questions/31051488/combining-two-matplotlib-colormaps colors_undersea = plt.cm.terrain(np.linspace(0, 0.17, 56)) colors_land = plt.cm.terrain(np.linspace(0.25, 1, 200)) ## combine them and build a new colormap colors = np.vstack((colors_undersea, colors_land)) cut_terrain_map = matplotlib.colors.LinearSegmentedColormap.from_list('cut_terrain', colors) plt.close('all') data = pd.read_csv('D:/Cumtb_Code/LeeCode/test.csv',names=['lon', 'lat','elev']) Long = data['lon']; Lat = data['lat']; Elev = data['elev'] Long = Long.values.tolist() Lat=Lat.values.tolist() Elev =Elev.values.tolist() tl = 0.5 tw = 0.5 lw = 1.5 S = 30 pts=1000000 [x,y] =np.meshgrid(np.arange(np.min(Long),np.max(Long),0.01),np.arange(np.min(Lat),np.max(Lat),0.01)) z = griddata((Long, Lat), Elev, (x, y), method='linear') x = np.matrix.flatten(x) y = np.matrix.flatten(y) z = np.matrix.flatten(z) cmap = plt.get_cmap('terrain') fig,ax = plt.subplots() norm = FixPointNormalize(sealevel=0,vmax=np.max(z)-400,vmin=np.min(z)+250) plt.scatter(x,y,1,z,cmap =cut_terrain_map,norm=norm) cbar = plt.colorbar(label='Elevation above sea level [m]') cbar.ax.tick_params(size=3,width =1) cbar.ax.tick_params(which='major',direction='in',bottom=True, top=True, left=False, right=True,length=tl*2.5,width=tw+1.5,color='k') cbar.outline.set_linewidth(lw) plt.xlabel('Longitude [°]') plt.ylabel('Latitude [°]') plt.gca().set_aspect('equal') plt.xlim(116,121) plt.ylim(30,35) # plt.xlim(40,52) # plt.ylim(-10,-28) plt.gca().invert_yaxis() plt.yticks([30,31,32,33,34,35]) # plt.yticks([-10,-15,-20,-25]) ax.set_yticks(np.linspace(30,35,6), minor=True) # ax.set_yticks(np.linspace(-28,-10,19), minor=True) ax.tick_params(colors='k') ax.spines['top'].set_linewidth(lw) ax.spines['bottom'].set_linewidth(lw) ax.spines['right'].set_linewidth(lw) ax.spines['left'].set_linewidth(lw) ax.tick_params(which='major',direction='in',bottom=True, top=True, left=True, right=True,length=tl*2,width=tw+1.5,color='k') ax.minorticks_on() ax.tick_params(which='minor',direction='in',bottom=True, top=True, left=True, right=True,color='k',length=tl+1,width=tw) plt.rcParams["font.family"] = "charter" plt.rcParams.update({'font.size': S}) ax.set_yticks(np.linspace(30,35,6), minor=True) # ax.set_yticks(np.linspace(-28,-10,19), minor=True) #显示 plt.savefig('D:/Cumtb_Code/LeeCode/plot9.png',dpi=800) plt.show()
输出结果
提示代码
''' Description: henggao_note version: v1.0.0 Date: 2022-04-01 10:41:56 LastEditors: henggao LastEditTime: 2022-04-02 08:58:30 ''' import pandas as pd import numpy as np from scipy.interpolate import griddata a = [1,2,3] b = [1,2,3] ans = [4,5,6,3,4,5,2,3,4] A,B = np.meshgrid(a,b) # print(A) # print(B) X_star = np.hstack((A.flatten()[:,None], B.flatten()[:,None])) # print(X_star) m = [1.5,2.5] n = [1.5,2.5] M,N = np.meshgrid(m,n) # print(M,N) U = griddata(X_star, ans, (M, N), method='cubic') print(U) # CSV_FILE_PATH = 'D:/Cumtb_Code/LeeCode/test.csv' # data = pd.read_csv(CSV_FILE_PATH,names=['lon', 'lat','elev']) # # print(data) # Long = data['lon']; Lat = data['lat']; Elev = data['elev'] # a = Long.values.tolist() # b=Lat.values.tolist() # ans =Elev.values.tolist() # # print(a) # # print(b) # print(len(ans)) # # A,B = np.meshgrid(a,b) # # X_star = np.hstack((A.flatten()[:,None], B.flatten()[:,None])) # # print(X_star) # # print(a,b) # [m,n] =np.meshgrid(np.arange(np.min(a),np.max(a),0.01),np.arange(np.min(b),np.max(b),0.01)) # M,N = np.meshgrid(m,n) # U = griddata((a,b), ans, (M, N), method='linear') # print(U)
拓展
将这个代码改进，绘制剖面图。
通过插值可以使得图形平滑过渡。
