GMT、Python、matlab、Arcgis在地图背景上添加散点的实现

在论文出图时，经常需要绘制地理区划图，而我们需要往其中添加点位坐标。下面我将分别介绍在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文件

绘制的效果

【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)

实现效果如下：

【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

【4】arcgis

如果是加载txt，或者csv数据，先显示xy坐标，然后加载背景图（GIA）

大家有什么问题，欢迎与我交流！