python气象信息系统工程 第八章

/小书匠/python

第8章 Python绘图基础

8.1 Matplotlib与cartopy基础知识

8.1.1 绘图结构

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
# 引入绘图库
import numpy as np
import matplotlib.pyplot as plt
# 绘图数据计算函数
def f(t):
    return np.exp(-t) * np.cos(2*np.pi*t)
# 生成两个子图对应的横轴数据
t1 = np.arange(0.0, 5.0, 0.1)
t2 = np.arange(0.0, 5.0, 0.02)
# 生成figure
plt.figure(1)
# 生成子图1(位置为两行一列排布的从上往下数第一个)
plt.subplot(211)
plt.plot(t1, f(t1), 'bo', t2, f(t2), 'k')
# 生成子图2(位置为两行一列排布的从上往下数第二个)
plt.subplot(212)
plt.plot(t2, np.cos(2*np.pi*t2), 'r--')
# 进行绘图
plt.show()

8.1.2 Figure、Axes与GeoAxes

1
2
3
4
5
6
7
8
9
10
11
12
13
14
import numpy as np
import matplotlib.pyplot as plt
fig=plt.figure(figsize=(12,6))
for i in range(6):
    ax=fig.add_subplot(2,3,i+1)
    ax.text(0.5,0.5,f'{i+1}', fontsize=20,ha='center',va='center')
    ax.text(0.5,0.2,f'ax=fig.add_subplot(2,3,{i+1})', fontsize=10,ha='center',va='center')
ax = fig.add_axes([0.3, 0.35, 0.3, 0.2])
ax.text(0.5,0.5,'ax=fig.add_axes([0.3, 0.35, 0.2, 0.2])', fontsize=10, ha='center', va='center')
# 这里的[0.3, 0.35, 0.2, 0.2]指Axes的左轴距离左边界0.3个Figure的宽度，
# 下轴距离下边界0.35个Figure的高度，Axes宽度为0.2个Figure的宽度，
# Axes的高度为0.2个Figure的高度
fig.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
f = xr.open_dataset('/home/mw/input/pythonbook9857/1980060106.nc')
z = f['z'].loc[:,500,:,:][0] #取出第一个时刻的500hPa的位势场
lat = f['latitude']
lon = f['longitude']
plt.figure(figsize=(5, 5))
#ax1为Axes
ax1 = plt.subplot(1,2,1)
c1 = ax1.contour(lon, lat, z)
ax1.set_title('Axes')
#ax2为GeoAxes
ax2 = plt.subplot(1,2,2, projection=ccrs.PlateCarree())
c2 = ax2.contour(lon, lat, z,transform=ccrs.PlateCarree())
ax2.set_title('GeoAxes')

Warning: ecCodes 2.21.0 or higher is recommended. You are running version 2.14.1





Text(0.5, 1.0, 'GeoAxes')

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
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker

f = xr.open_dataset('/home/mw/input/pythonbook9857/1980060106.nc')
z = f['z'].loc[:,500,:,:][0]
lat = f['latitude']
lon = f['longitude']
plt.figure(figsize=(8, 6))
ax1 = plt.subplot(1,1,1, projection=ccrs.PlateCarree())
ax1.set_extent([60,180,0,90]) #设置GeoAxes的范围，4个数字分别代表左、右、上、下边界的经纬度
ax1.gridlines() #添加栅格线
#ax1.coastlines() #添加海岸线
#ax1.add_feature(cfeature.LAND) #添加大陆特征
#添加刻度
ax1.set_xticks(np.arange(60,210,30), crs=ccrs.PlateCarree())
ax1.set_yticks(np.arange(0,120,30), crs=ccrs.PlateCarree())
ax1.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax1.yaxis.set_major_formatter(cticker.LatitudeFormatter())
#绘制等值线
c1 = ax1.contour(lon, lat, z,transform=ccrs.PlateCarree())
plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
fig = plt.figure(figsize=[10, 5])
ax1 = fig.add_subplot(1, 2, 1,projection=ccrs.PlateCarree())
ax1.set_xticks(np.arange(-180,240,60), crs=ccrs.PlateCarree())
ax1.set_yticks(np.arange(-90,120,30), crs=ccrs.PlateCarree())
ax1.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax1.yaxis.set_major_formatter(cticker.LatitudeFormatter())
ax1.gridlines() #添加栅格线
# ax1.add_feature(cfeature.COASTLINE.with_scale('110m'))
ax2 = fig.add_subplot(1, 2, 2, projection=ccrs.PlateCarree(central_longitude = 120))
ax2.set_xticks(np.arange(-180,240,60), crs=ccrs.PlateCarree())
ax2.set_yticks(np.arange(-90,120,30), crs=ccrs.PlateCarree())
ax2.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax2.yaxis.set_major_formatter(cticker.LatitudeFormatter())
ax2.gridlines() #添加栅格线
# ax2.add_feature(cfeature.COASTLINE.with_scale('110m'))
plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
fig = plt.figure(figsize=[10, 5])
ax1 = fig.add_subplot(1, 2, 1, projection=ccrs.LambertConformal())
gl1 = ax1.gridlines(draw_labels=True,x_inline=False, y_inline=False) #添加栅格线，绘制刻度标签，但禁止在栅格线内绘制标签
gl1.rotate_labels = False #禁止刻度标签旋转
# ax1.add_feature(cfeature.COASTLINE.with_scale('110m'))
ax2 = fig.add_subplot(1, 2, 2, projection=ccrs.LambertConformal(central_longitude=120,cutoff=0))
gl2 = ax2.gridlines(draw_labels=True,x_inline=False, y_inline=False) #添加栅格线
gl2.rotate_labels = False
# ax2.add_feature(cfeature.COASTLINE.with_scale('110m'))
plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
fig = plt.figure(figsize=[10, 5])
ax1 = fig.add_subplot(1, 2, 1, projection=ccrs.Mercator())
# ax1.add_feature(cfeature.COASTLINE.with_scale('110m'))
gl1 = ax1.gridlines(draw_labels=True,x_inline=False, y_inline=False) #添加栅格线
#隐藏上、右坐标轴刻度
gl1.xlabels_top = False
gl1.ylabels_right = False
plt.show()

/opt/conda/lib/python3.7/site-packages/cartopy/mpl/gridliner.py:307: UserWarning: The .xlabels_top attribute is deprecated. Please use .top_labels to toggle visibility instead.
  warnings.warn('The .xlabels_top attribute is deprecated. Please '
/opt/conda/lib/python3.7/site-packages/cartopy/mpl/gridliner.py:343: UserWarning: The .ylabels_right attribute is deprecated. Please use .right_labels to toggle visibility instead.
  warnings.warn('The .ylabels_right attribute is deprecated. Please '

1
2
3
4
5
6
7
8
9
10
11
12
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
fig = plt.figure(figsize=[10, 5])
ax1 = fig.add_subplot(1, 2, 1, projection=ccrs.NorthPolarStereo())
# ax1.add_feature(cfeature.COASTLINE.with_scale('110m'))
gl1 = ax1.gridlines(draw_labels=True,x_inline=False, y_inline=True) #添加栅格线
ax2 = fig.add_subplot(1, 2, 2, projection=ccrs.SouthPolarStereo())
# ax2.add_feature(cfeature.COASTLINE.with_scale('110m'))
gl2 = ax2.gridlines(draw_labels=True,x_inline=False, y_inline=True) #添加栅格线
plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
import matplotlib.path as mpath
import matplotlib.pyplot as plt
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
fig = plt.figure(figsize=[5, 5])
ax1 = fig.add_subplot(1, 1, 1, projection=ccrs.NorthPolarStereo())
ax1.set_extent([-180, 180, 30, 90], ccrs.PlateCarree())
gl1 = ax1.gridlines(draw_labels=True,x_inline=False, y_inline=False) #添加栅格线
gl1.rotate_labels = False #禁止标签旋转，设置为True时标签与栅格线平行
# ax1.add_feature(cfeature.COASTLINE.with_scale('50m'))
plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
import matplotlib.path as mpath
import matplotlib.pyplot as plt
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
fig = plt.figure(figsize=[5, 5])
ax1 = fig.add_subplot(1, 1, 1, projection=ccrs.NorthPolarStereo())
ax1.set_extent([-180, 180, 30, 90], ccrs.PlateCarree())
ax1.gridlines()
ax1.add_feature(cfeature.COASTLINE.with_scale('50m'))
# 生成一个圆形的Path
theta = np.linspace(0, 2*np.pi, 100)
center, radius = [0.5, 0.5], 0.5
verts = np.vstack([np.sin(theta), np.cos(theta)]).T
circle = mpath.Path(verts * radius + center)
# 将该Path设置为GeoAxes的边界
ax1.set_boundary(circle, transform=ax1.transAxes)
plt.show()

/opt/conda/lib/python3.7/site-packages/cartopy/io/__init__.py:260: DownloadWarning: Downloading: https://naciscdn.org/naturalearth/50m/physical/ne_50m_coastline.zip
  warnings.warn('Downloading: {}'.format(url), DownloadWarning)

8.2 地理绘图基础

8.2.2 在GeoAxes上绘制

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.mpl.ticker as cticker
import geopandas as gpd

import geopandas as gpd
gdf = gpd.GeoDataFrame.from_file('/home/mw/input/pythonbook9857/circle.json', encoding='utf8')

fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
ax.set_extent([115.5, 123, 29.5, 36.5], ccrs.PlateCarree())
ax.add_geometries(gdf['geometry'], crs=ccrs.PlateCarree(), facecolor='none', edgecolor='black')
# 设置轴刻度
ax.set_xticks(np.arange(116, 123, 1), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(30, 36, 1), crs=ccrs.PlateCarree())
ax.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax.yaxis.set_major_formatter(cticker.LatitudeFormatter())
plt.show()

8.2.3 几何数据筛选示例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.mpl.ticker as cticker
import geopandas as gpd

name_a = gdf.loc[gdf['name']=='A']

fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
ax.set_extent([115.5, 123, 29.5, 36.5], ccrs.PlateCarree())
ax.add_geometries(name_a['geometry'], crs=ccrs.PlateCarree(), facecolor='none', edgecolor='black')
# 设置轴刻度
ax.set_xticks(np.arange(116, 123, 1), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(30, 36, 1), crs=ccrs.PlateCarree())
ax.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax.yaxis.set_major_formatter(cticker.LatitudeFormatter())
plt.show()

8.2.4 多边形合并

1
2
3
4
5
6
7
8
9
10
11
12
import matplotlib.pyplot as plt
from shapely.ops import unary_union
import geopandas as gpd
import cartopy.crs as ccrs
# 读取GeoJSON文件
gdf = gpd.GeoDataFrame.from_file('/home/mw/input/pythonbook9857/circle.json', encoding='utf8') #也可以是shapefile
# 合并多边形
geom = unary_union([geom if geom.is_valid else geom.buffer(0) for geom in gdf['geometry']])

print(type(gdf['geometry']), len(gdf['geometry']))
print(type(geom))

<class 'geopandas.geoseries.GeoSeries'> 2
<class 'shapely.geometry.polygon.Polygon'>

8.3 颜色表（colormap）

8.3.3 创建自定义色标

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
# 样例数据生成，此段可以忽略
N = 100
X, Y = np.mgrid[-3:3:complex(0, N), -2:2:complex(0, N)]
Z1 = np.exp(-X**2 - Y**2)
Z2 = np.exp(-(X - 1)**2 - (Y - 1)**2)
Z = (Z1 - Z2) * 2
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
levels = [0.1, 0.3, 0.5, 0.9]
cmap = mcolors.ListedColormap(['r', 'g', 'b'])
cmap.set_under('c')
cmap.set_over('m')
norm = mcolors.BoundaryNorm(levels, cmap.N)
pcm = ax.contourf(Z, levels=levels, cmap=cmap, norm=norm, extend='both')
fig.colorbar(pcm, ax=ax, extend='both')
plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
N = 100
X, Y = np.mgrid[-3:3:complex(0, N), -2:2:complex(0, N)]
Z1 = np.exp(-X**2 - Y**2)
Z2 = np.exp(-(X - 1)**2 - (Y - 1)**2)
Z = (Z1 - Z2) * 2
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
levels = [0.1, 0.3, 0.5, 0.9]
cmap, norm = mcolors.from_levels_and_colors(levels, ['c', 'r', 'g', 'b', 'm'], 'both')
pcm = ax.contourf(Z, levels=levels, cmap=cmap, norm=norm, extend='both')
fig.colorbar(pcm, ax=ax, extend='both')
plt.show()

8.4 图像显示与保存

8.4.1 图像显示

1
2
3
4
5
6
7
8
9
10
11
12
import numpy as np
import matplotlib.pyplot as plt
#生成正弦函数
x = np.linspace(0, 10, 500)
y = np.sin(x)
#绘制图像
fig = plt.figure(figsize=(10, 8))
ax1 = fig.add_subplot(1,1,1)
ax1.plot(x, y)
#显示图像
plt.show()

8.4.2 图像保存

1
2
3
4
5
6
7
8
9
10
import numpy as np
import matplotlib.pyplot as plt
#生成正弦函数
x = np.linspace(0, 10, 500)
y = np.sin(x)
#绘制图像
fig = plt.figure(figsize=(10, 8))
ax1 = fig.add_subplot(1,1,1)
ax1.plot(x, y)
plt.savefig('./sinx.png', dpi=400, bbox_inches='tight')