产品
在论文出图时，经常需要绘制地理区划图，而我们需要往其中添加点位坐标。下面我将分别介绍在GMT、Python、matlab、Arcgis中添加散点的方法。以绘制全球的GPS站点为例。
【1】GMT
实现的代码如下：
gmtset MAP_FRAME_TYPE plain
set R=-180/180/-90/90
set J=Q25c
set PS=IGS_global.ps
REM gmt xyz2grd b.txt -R%R% -I1 -Gtmp.grd
REM gmt grdsample tmp.grd -Gtmp.grd -I0.01
REM gmt makecpt -Cseis -T10/60/0.1 -D >tem.cpt
REM gmt psclip zhujiang.txt -J%J% -R%R% -Am -K >%PS%
REM gmt grdimage tmp.grd -R%R% -J%J% -Ctem.cpt -B1f1/1f1WenS -Xc -Yc -K>%PS%
REM gmt grdgradient ETOPO1_Bed_g_gdal.grd -Ne0.7 -A50 -Gtmp2.grd
gmt makecpt -Cetopo1 -D >tem1.cpt
gmt grdimage ETOPO1_Bed_g_gdal.grd -Itmp2.grd -R%R% -J%J% -Ctem1.cpt -K -Xc -Yc>%PS%
gmt pscoast -R%R% -J%J% -A1000 -N1 -I1 -Wthinnest -BWeSn -Xc -Yc -O -K>> %ps%
gmt psxy ddd.txt -R%R% -J%J% -Sc0.04c -W0.01,red -Gred -O -K >> %ps%
gmt psbasemap -R%R% -J%J% -Bf30a30 -BWeSn -Bg60 -Xc -Yc --FONT_ANNOT_PRIMARY=16p,Helvetica,black -O >>%PS%
gmt psconvert %PS% -A -Tg -P
其中的ddd.txt文件
绘制的效果
GMT
【2】Python
import os import numpy as np import matplotlib as mpl import matplotlib.pyplot as plt import cartopy.crs as ccrs from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER from mpl_toolkits.basemap import Basemap import pandas as pd import h5py FILE_NAME = './dataset/ATL08_20210601183700_10531105_004_01.h5' with h5py.File(FILE_NAME, mode='r') as f: latvar = f['/gt1l/land_segments/latitude'] lat = latvar[:] lonvar = f['/gt1l/land_segments/longitude'] lon = lonvar[:] dset_name = '/gt1l/land_segments/dem_h' datavar = f[dset_name] data = datavar[:] units = datavar.attrs['units'] long_name = datavar.attrs['long_name'] _FillValue = datavar.attrs['_FillValue'] # Handle FillValue data[data == _FillValue] = np.nan data = np.ma.masked_where(np.isnan(data), data) # Find the middle location. lat_m = lat[int(lat.shape[0]/2)] lon_m = lon[int(lon.shape[0]/2)] # print(lon_m) # Let's use ortho projection. # orth = ccrs.Orthographic(central_longitude=lon_m, # central_latitude=lat_m, # globe=None) # ax = plt.axes(projection=orth) m = Basemap(projection='cyl',lat_0=lat_m, lon_0=lon_m) # Remove the following line to see a zoom-in view. # ax.set_global() # Put grids. # gl = ax.gridlines(draw_labels=True, dms=True) m.drawparallels(np.arange(-90., 90., 20.), labels=[1,0,0,0], fontsize=10) m.drawmeridians(np.arange(-180., 180., 30.), labels=[0,0,0,1], fontsize=10) m.drawcoastlines() m.shadedrelief(scale=0.1) #添加山体阴影的全球地图 # Put coast lines. # ax.coastlines() # m.drawcoastlines() # Plot on map. # n = 51 # x = np.random.randn(n) * 180 # print(x) # y = np.random.randn(n) * 90 # print(y) data = pd.read_table('./dataset/ddd.txt') a=np.loadtxt('./dataset/ddd.txt') print(a) m.scatter(a[:,0], a[:,1], marker='v', s=100, facecolor='#00BFFF', edgecolor='k', linewidth=1) m.scatter(lon,lat,marker = 'o',s = 10,facecolor = 'red') # p = plt.scatter(lon, lat, c=data, s=1, cmap=plt.cm.jet, # transform=ccrs.PlateCarree()) # Put grid labels at left and bottom only. # gl.top_labels = False # gl.right_labels = False # Put degree N/E label. # gl.xformatter = LONGITUDE_FORMATTER # gl.yformatter = LATITUDE_FORMATTER # Adjust colorbar size and location using fraction and pad. # cb = plt.colorbar(p, fraction=0.022, pad=0.01) units = units.decode('ascii', 'replace') # cb.set_label(units, fontsize=8) basename = os.path.basename(FILE_NAME) long_name = long_name.decode('ascii', 'replace') plt.title('{0}\n{1}'.format(basename, long_name)) plt.show() fig = plt.gcf() pngfile = "{0}.py.png".format(basename) fig.savefig(pngfile)
实现效果如下：
python
【3】matlab
实现代码如下：
[ grid_pgr ] = gmt_readpgr( 'GIA_n100_mass_0km.txt' );
% get delta-SA data from the TEOS-10 gsw atlas at 2500 dbar
[LG,LT]=meshgrid(0.5:360,-89.5:90);
dSA=ones(size(LG));
grid_pgr= flipud(grid_pgr);
dSA(:)=grid_pgr;
% Rearrange data to lie in the longitude limits I give for the
% % projection
% ind=[31:361 1:30]; % Move left side to right
% dSA=dSA(:,ind);
% LT=LT(:,ind);
% LG=LG(:,ind);LG(LG>30)=LG(LG>30)-360; %...and subtract 360 to some longitudes
% clf;
m_proj('robinson','lon',[0 360]);
m_pcolor(LG,LT,dSA);
m_grid('tickdir','out','linewi',2);
m_coast('color','k');
A=textread('ddd.txt');
for i=1:length(A)
[X,Y]=m_ll2xy(A(i,1),A(i,2));%convert cordinates
line(X,Y,'marker','.','markersize',8,'color','r');
end
[X,Y]=m_ll2xy(106,30);
% This is a perceptually uniform jet-like color scale, but in m_colmap
% we can add some simple graduated steps to make the pcolor look a little
% more like a contourf
% colormap(m_colmap('land','step',10));
% colormap(m_colmap('blue','step',10));
% colormap(m_colmap('odv','step',10));
% colormap(m_colmap('gland','step',10));
colormap(m_colmap('diverging','step',5));
% m_northarrow(-123 -4.5/60,49+19.5/60,1/60,'type',4,'aspect',1.5);
% hold on;
% [cs,h]=contour(LG,LT,dSA,'color','k');
% clabel(cs,h,'fontsize',6);
% h=colorbar('northoutside');
% title(h,'GIA global','fontsize',14);
% set(h,'pos',get(h,'pos')+[.2 .05 -.4 0],'tickdir','out')
set(gcf,'color','w'); % Need to do this otherwise 'print' turns the lakes black
matlab
【4】arcgis
如果是加载txt，或者csv数据，先显示xy坐标，然后加载背景图（GIA）
ArcGIS
大家有什么问题，欢迎与我交流！
155-2731-8020
