http://www.gnss.help/2018/04/24/cartopy-gallery/
具体绘制的时候,还可以参考https://scitools.org.uk/cartopy/docs/latest/gallery/index.html
绘制中国政区图
本图使用的中国国界、省界、十段线数据中,边界线数据块使用 “>” 号分隔。因此首先将其内容按照 “>” 切块,然后加载到 NumPy 中。绘图使用的代码如下:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 |
import numpy as np import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature # Load the border data, CN-border-La.dat is downloaded from # http://gmt-china.org/datas/CN-border-La.dat with open('CN-border-La.dat') as datf: context = datf.read() blocks = [cnt for cnt in context.split('>') if len(cnt) > 0] borders = [np.fromstring(block, dtype=float, sep=' ') for block in blocks] # Set figure size fig = plt.figure(figsize=[10, 8]) # Set projection and plot the main figure ax = plt.axes(projection=ccrs.LambertConformal(central_latitude=90, central_longitude=105)) # Add ocean, land, rivers and lakes ax.add_feature(cfeature.OCEAN.with_scale('50m')) ax.add_feature(cfeature.LAND.with_scale('50m')) ax.add_feature(cfeature.RIVERS.with_scale('50m')) ax.add_feature(cfeature.LAKES.with_scale('50m')) # Plot border lines for line in borders: ax.plot(line[0::2], line[1::2], '-', color='gray', transform=ccrs.Geodetic()) # Plot gridlines ax.gridlines(linestyle='--') # Set figure extent ax.set_extent([80, 130, 13, 55]) # Plot South China Sea as a subfigure sub_ax = fig.add_axes([0.741, 0.11, 0.14, 0.155], projection=ccrs.LambertConformal(central_latitude=90, central_longitude=115)) # Add ocean, land, rivers and lakes sub_ax.add_feature(cfeature.OCEAN.with_scale('50m')) sub_ax.add_feature(cfeature.LAND.with_scale('50m')) sub_ax.add_feature(cfeature.RIVERS.with_scale('50m')) sub_ax.add_feature(cfeature.LAKES.with_scale('50m')) # Plot border lines for line in borders: sub_ax.plot(line[0::2], line[1::2], '-', color='gray', transform=ccrs.Geodetic()) # Set figure extent sub_ax.set_extent([105, 125, 0, 25]) # Show figure plt.show() |
绘制 IGS 站点分布图
本图使用的 IGS 核心站与 MGEX 项目站点,及其坐标均来自 IGS 网站。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 |
import numpy as np import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature # Load the coordinate of IGS Core & MGEX sites, The CSV files are # exported from: http://www.igs.org/network # Download igs-core.csv & mgex.csv from # http://ocijpaoi3.bkt.clouddn.com/igs-core.csv # http://ocijpaoi3.bkt.clouddn.com/mgex.csv igs_core = np.recfromcsv('igs-core.csv', names=True, encoding='utf-8') mgex = np.recfromcsv('mgex.csv', names=True, encoding='utf-8') fig = plt.figure(figsize=[9, 6]) # Set projection ax = plt.axes(projection=ccrs.Robinson()) # Add ocean and land ax.add_feature(cfeature.LAND) ax.add_feature(cfeature.OCEAN) # Add MGEX & IGS core sites ax.plot(mgex['longitude'], mgex['latitude'], 'o', color='tomato', label='MGEX', transform=ccrs.Geodetic()) ax.plot(igs_core['longitude'], igs_core['latitude'], '*', color='darkmagenta', label='IGS Core', transform=ccrs.Geodetic()) # Plot gridlines ax.gridlines(linestyle='--') # Set figure extent ax.set_global() # Set legend location plt.legend(loc='lower right') # Show figure plt.show() |
绘制 GNSS 控制网
这里使用的 IGS 站点坐标数据同样来自 IGS 网站。我已将其整理成一个 CSV 格式的文件:euro-igs,你可以直接下载使用。这里使用 matplotlib.tri
中的 Triangulation
来根据输入的点位坐标来创建 Delaunay 三角网,然后使用 plt.triplot()
方法绘制。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 |
import numpy as np import matplotlib.pyplot as plt import matplotlib.tri as tri import cartopy.crs as ccrs import cartopy.feature as cfeature # Load coordinate of the IGS sites in Europe, this CSV file is # exported from: http://www.igs.org/network # Download euro-igs.csv from: # http://ocijpaoi3.bkt.clouddn.com/euro-igs.csv igs_sites = np.recfromcsv('euro-igs.csv', names=True, encoding='utf-8') # Generate Delaunay triangles triangles = tri.Triangulation(igs_sites['longitude'], igs_sites['latitude']) fig = plt.figure(figsize=[6, 8]) # Set projection ax = plt.axes(projection=ccrs.LambertConformal(central_latitude=90, central_longitude=10)) # Add ocean, land, rivers and lakes ax.add_feature(cfeature.OCEAN.with_scale('50m')) ax.add_feature(cfeature.LAND.with_scale('50m')) ax.add_feature(cfeature.RIVERS.with_scale('50m')) ax.add_feature(cfeature.LAKES.with_scale('50m')) # Plot triangles plt.triplot(triangles, transform=ccrs.Geodetic(), marker='o', linestyle='-') # Plot gridlines ax.gridlines(linestyle='--') # Set figure extent ax.set_extent([-10, 30, 30, 73]) # Show figure plt.show() |
突出显示某区域
Cartopy 默认对所有地物使用相同的外观,如果需要突出显示某些地物,就必须进行筛选。这里从 Natural Earth 提供的小比例尺国界数据中,提取出欧盟国家,然后使用 ax.add_geometries()
方法将它们加入到绘图元素中。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 |
import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature import cartopy.io.shapereader as shpreader # Country names in Europe Union, exported from Wikipedia: # https://en.wikipedia.org/wiki/European_Union eu_country_names = ('Austria', 'Belgium', 'Bulgaria', 'Croatia', 'Cyprus', 'Czechia', 'Denmark', 'Estonia', 'Finland', 'France', 'Germany', 'Greece', 'Hungary', 'Ireland', 'Italy', 'Latvia', 'Lithuania', 'Luxembourg', 'Malta', 'Netherlands', 'Poland', 'Portugal', 'Romania', 'Slovakia', 'Slovenia', 'Spain', 'Sweden', 'United Kingdom') # Create a .shp file reader shp_file = shpreader.natural_earth(resolution='110m', category='cultural', name='admin_0_countries') shp_reader = shpreader.Reader(shp_file) # Reader the shape file records = shp_reader.records() # Collect the European Union countries eu_countries = [] for rec in records: if rec.attributes['NAME'] in eu_country_names: eu_countries.append(rec) # Start ploting fig = plt.figure(figsize=[6, 6]) # Set projection ax = plt.axes(projection=ccrs.Orthographic(central_latitude=50, central_longitude=10)) # Add land ax.add_feature(cfeature.LAND, facecolor='lightgray') # Plot the European Union countries for country in eu_countries: ax.add_geometries(country.geometry, crs=ccrs.PlateCarree(), facecolor='darkgreen') # Plot gridlines ax.gridlines(linestyle='--') # Show figure plt.show() |