在論文出圖時(shí)，經(jīng)常需要繪制地理區(qū)劃圖，而我們需要往其中添加點(diǎn)位坐標(biāo)。下面我將分別介紹在GMT、Python、matlab、Arcgis中添加散點(diǎn)的方法。以繪制全球的GPS站點(diǎn)為例。

【1】GMT

實(shí)現(xiàn)的代碼如下：

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)

實(shí)現(xiàn)效果如下：

【3】matlab

實(shí)現(xiàn)代碼如下：

[ 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數(shù)據(jù)，先顯示xy坐標(biāo)，然后加載背景圖（GIA）

大家有什么問(wèn)題，歡迎與我交流！