Basemap系列教程：绘图

bugsuse

发布于 2020-04-21 17:46:00

4K0

发布于 2020-04-21 17:46:00

文章被收录于专栏：气象杂货铺气象杂货铺

说了这么久，也看了这么多例子了。这一次我们看一下一些绘图命令。

annotate

创建含有箭头的文本以表示感兴趣的点。要创建无箭头的文本，可使用text。

annotate(*args, **kwargs)

text 方法并不属于 Basemap，而是属于 matplotlib。因此，必须要从 plot 或 axis 实例进行调用
第一个参数是文本字符串
xy 列表是箭头所指向的点的 x 和 y 坐标。依赖于 xycoords 参数设置。
xycoords 表示 xy 使用的坐标类型 1） data 意味着使用的坐标和 data 使用的坐标一致 2） offset point 意味着点的偏移 3） axes pixels 表示从轴左下角开始的像素 4）其余选项可在 annotate 文档查看
xytext 表示文本的 x 和 y 坐标，同时也是箭头的起始点
textcoords 是 xytext所使用的坐标轴类型，和 xycoords 具有相同的选项
arrowprops 表示箭头的属性，为可选参数。详细解释见 Line2D
color 表示文本的颜色

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt


map = Basemap(projection='ortho', lat_0=0, lon_0=0)

map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='coral',lake_color='aqua')
map.drawcoastlines()


x, y = map(2, 41)
x2, y2 = (-90, 10)

plt.annotate('Barcelona', xy=(x, y),  xycoords='data', xytext=(x2, y2), textcoords='offset points', color='r', arrowprops=dict(arrowstyle="fancy", color='g'))

x2, y2 = map(0, 0)
plt.annotate('Barcelona', xy=(x, y),  xycoords='data', xytext=(x2, y2), textcoords='data', arrowprops=dict(arrowstyle="->"))
plt.show()

barbs

在地图中绘制风杆

barbs(x, y, u, v, *args, **kwargs)

关于 barb 的详细说明可以查看 matplotlib 文档。[注1]

x 和 y 是给定的网格点的位置，如果 latlon 参数为 True，这些值将被假设为地理学坐标中的点，否则视为地图坐标系中的点
u 和 v 是以 knot 为单位的左右和上下方向风的大小。如果投影为非圆柱投影（即除了 cyl，merc，gall，mill [注2]），应该使用 rotate_vector 或 transform_scalar 对 u 和 v 进行旋转

同时还有很多其它可选参数 [注3]。一些常用参数如下：

pivot ：改变 barbs 的旋转点，默认为 'tip'，可设置为 'middle'
barbcolor ：可以改变颜色。如果使用序列，颜色将被交替设置
fill_empty：如果风速小于 5 kt 时，填充圆
barb_increments ：用于改变 barb 添加的 tick 和 flag 的值

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from osgeo import gdal
import numpy as np


map = Basemap(llcrnrlon=-93.7, llcrnrlat=28., urcrnrlon=-66.1, urcrnrlat=39.5,
              projection='lcc', lat_1=30., lat_2=60., lat_0=34.83158, lon_0=-98.)

ds = gdal.Open("../sample_files/wrf.tiff")
lons = ds.GetRasterBand(4).ReadAsArray()
lats = ds.GetRasterBand(5).ReadAsArray()
u10 = ds.GetRasterBand(1).ReadAsArray()
v10 = ds.GetRasterBand(2).ReadAsArray()

x, y = map(lons, lats)

yy = np.arange(0, y.shape[0], 4)
xx = np.arange(0, x.shape[1], 4)

points = np.meshgrid(yy, xx)

map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='#cc9955', lake_color='aqua', zorder = 0)
map.drawcoastlines(color = '0.15')

map.barbs(x[points], y[points], u10[points], v10[points], 
    pivot='middle', barbcolor='#333333')


plt.show()

使用 barbs时最主要的问题是点的密度可能会很大，而这一方法又不能自动忽略这些点 1）可仅使用矩阵数据的 1/4 进行绘图 2）矩阵中的点要含有需要的数据点 3）传递给 barbs 的参数可使用点阵进行选取，因此可每 4 个元素选择 1 个值
u 和 v 被认为是地图投影坐标，而不是经纬度坐标。如果是地理学坐标系的话，可使用 rotate_vector 方法进行适当的旋转
pivot 参数可从 middle 进行 barbs 旋转，因为当风旋转时，从上部旋转（pivot = ‘tip’）的效果会产生奇怪的现象

contour

创建等值线图

contour(x, y, data)

x 和 y 是和 data 矩阵具有相同大小的矩阵，包含地图坐标系中元素的位置
data 是包含要绘图的数据的矩阵
第四个参数可以被传递，包含绘制等值线图时的一系列等值线值
默认的 colormap 是 jet，但可通过设置 cmap 参数改变 colormap
参数 tri = True 时，网格会被认为是不规则网格，效果差异可在此链接进行查看 [注4]
其余可能的参数可查看matplotlib 官方文档 [注5]
labels 可以被添加到图中，比如 Basemap系列教程：基本函数 部分 contour 例子说明

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from osgeo import gdal
from numpy import linspace
from numpy import meshgrid

map = Basemap(projection='tmerc', 
              lat_0=0, lon_0=3,
              llcrnrlon=1.819757266426611, 
              llcrnrlat=41.583851612359275, 
              urcrnrlon=1.841589961763497, 
              urcrnrlat=41.598674173123)

ds = gdal.Open("../sample_files/dem.tiff")
data = ds.ReadAsArray()

x = linspace(0, map.urcrnrx, data.shape[1])
y = linspace(0, map.urcrnry, data.shape[0])

xx, yy = meshgrid(x, y)

map.contour(xx, yy, data)

plt.show()

contourf

创建填充等值线图

contourf(x, y, data)

x 和 y 是和 data 矩阵具有相同大小的矩阵，包含地图坐标系中元素的位置
data 是包含要绘图的数据的矩阵
第四个参数可以被传递，包含绘制等值线图时的一系列等值线值
默认的 colormap 是 jet，但可通过设置 cmap 参数改变 colormap
参数 tri = True 时，网格会被认为是不规则网格，效果差异可在此链接进行查看 [注4]
其余可能的参数可查看matplotlib 官方文档 [注5]

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from osgeo import gdal
from numpy import linspace
from numpy import meshgrid

map = Basemap(projection='tmerc', 
              lat_0=0, lon_0=3,
              llcrnrlon=1.819757266426611, 
              llcrnrlat=41.583851612359275, 
              urcrnrlon=1.841589961763497, 
              urcrnrlat=41.598674173123)

ds = gdal.Open("../sample_files/dem.tiff")
data = ds.ReadAsArray()

x = linspace(0, map.urcrnrx, data.shape[1])
y = linspace(0, map.urcrnry, data.shape[0])

xx, yy = meshgrid(x, y)

map.contourf(xx, yy, data)

plt.show()

hexbin

由点集绘制六角箱。可以绘制每一个六边形出现的次数，或是给出每一次出现的权重。

hexbin(x, y, **kwargs)

更多的信息可以查看官方文档。[注6]

x 和 y 是包含点的坐标的numpy数组，常规的 list 不行，只能是 numpy 数组
gridsize 设置的是 x 方向六边形（bin）的数目。默认是 100
C 参数对应的是每一个点的值的数组。这些值默认由 numpy.mean 函数在每个 bin（六边形）进行处理
reduce_C_function 是对每一个 bin（六边形）的元素所执行的函数。默认是 numpy.mean() 函数
bins 函数可以控制计数函数的行为
mincnt 表示绘制每一个 bin（六边形）出现的最小值。默认为0，为了防止绘制六边形时没有出现，可设置为1
cmap 表示设置 colormap
edgecolors 表示六边形边缘的颜色。为了看到结果，linewidth 不能是 None
linewidths 六边形的线宽，默认为 None，即不绘制边界
alpha 表示 layer的透明度

注意：旧版本的库不支持 hexbin

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import matplotlib.colors as colors
from numpy import array
from numpy import max


map = Basemap(llcrnrlon=-0.5,llcrnrlat=39.8,urcrnrlon=4.,urcrnrlat=43.,
             resolution='i', projection='tmerc', lat_0 = 39.5, lon_0 = 1)


map.readshapefile('../sample_files/lightnings', 'lightnings')

x = []
y = []
c = []

for info, lightning in zip(map.lightnings_info, map.lightnings):
    x.append(lightning[0])
    y.append(lightning[1])
    
    if float(info['amplitude']) < 0:
        c.append(-1 * float(info['amplitude']))
    else:
        c.append(float(info['amplitude']))
    
plt.figure(0)

map.drawcoastlines()
map.readshapefile('../sample_files/comarques', 'comarques')

map.hexbin(array(x), array(y))

map.colorbar(location='bottom')



plt.figure(1)

map.drawcoastlines()
map.readshapefile('../sample_files/comarques', 'comarques')

map.hexbin(array(x), array(y), gridsize=20, mincnt=1, cmap='summer', bins='log')

map.colorbar(location='bottom', format='%.1f', label='log(# lightnings)')


plt.figure(2)

map.drawcoastlines()
map.readshapefile('../sample_files/comarques', 'comarques')

map.hexbin(array(x), array(y), gridsize=20, mincnt=1, cmap='summer', norm=colors.LogNorm())

cb = map.colorbar(location='bottom', format='%d', label='# lightnings')

cb.set_ticks([1, 5, 10, 15, 20, 25, 30])
cb.set_ticklabels([1, 5, 10, 15, 20, 25, 30])


plt.figure(3)

map.drawcoastlines()
map.readshapefile('../sample_files/comarques', 'comarques')

map.hexbin(array(x), array(y), C = array(c), reduce_C_function = max, gridsize=20, mincnt=1, cmap='YlOrBr', linewidths=0.5, edgecolors='k')

map.colorbar(location='bottom', label='Mean amplitude (kA)')


plt.show()

使用四个不同的例子来展示不同的选项
14 - 25 行读取所有的点，这些点是和闪电相关的数据。这些在 Basemap系列教程：使用shapefiles绘制地图读取点数据 部分进行了说明
第一个例子展示的是 hexbin 最小值的使用。注意： bin 非常小，那些出现次数为0的区域也被绘制了（深蓝色区域），而且shapefile 数据边界框外的区域（白色区域）也没有数据
第二个例子使用 gridsize 参数改变了 bin 的大小，使用 mincnt设置出现次数的最小值，从而进行绘图，并且设置 bins 为'log'，从而使颜色以对数尺度显示。注意：color bar 显示的是对数值
为了避免 colorbar 显示为对数值，可以使用 stackoverflow 中提到的方法[注7] 替换 bins 参数，但仍使用相同的数据绘图，但可以得到一个更好的colorbar。（译注：即图3所显示）
最后一个例子展示了如何使用 C 参数 1）闪电的振幅模块被存储在 c 变量中，然后传递给 C 参数

（1） reduce_C_function 参数用于显示每一个 bin 的最大值，从而代替平均值

（2） linewidths 和 edgecolors 使每一个六边形的边界被绘制

基本用法

对数刻度，使用不同的 hexagon 大小

对数刻度，具有更合适的colorbar

使用 C 参数，并且绘制六边形边界

imshow

在地图上绘制图像。图像可以是常规的 rgb 图，也可以是用 cmap 填充的图。

imshow(*args, **kwargs)

详细文档可查看 matplotlib 官方文档。[注8]

第一个参数是图像数组。如果有3个 band 的话，将被认为是 RGB图像，然后绘制图像。如果仅有 1 band，根据 cmap 参数（默认为 jet）将创建伪彩色。当数组具有2个或超过3个 band 时将会触发异常
extent 参数用于设置地图坐标中图像四个拐角的位置。必须是一个包含(x0, x1, y0, y1)的序列
origin 参数改变图像初始线的位置。默认，第一条线的位置由 image.origin 定义，可以使用 'upper' 或 ‘lower’ 改变此值
如果只有 1 个 band的话， cmap 用于设置 colormap
alpha 用于设置透明度，其值在 [0 1]之间
interpolation 参数用于当需要重置图像大小时改变插值模式。可能的值有‘none’, ‘nearest’, ‘bilinear’, ‘bicubic’, ‘spline16’, ‘spline36’, ‘hanning’, ‘hamming’, ‘hermite’, ‘kaiser’, ‘quadric’, ‘catrom’, ‘gaussian’, ‘bessel’, ‘mitchell’, ‘sinc’, ‘lanczos‘

第一个例子显示了使用 jpg 图像和一个域作为伪彩色图像：

from mpl_toolkits.basemap import Basemap
from mpl_toolkits.basemap import shiftgrid
import matplotlib.pyplot as plt
from osgeo import gdal
import numpy as np

map = Basemap(projection='tmerc', 
              lat_0=0, lon_0=3,
              llcrnrlon=1.819757266426611, 
              llcrnrlat=41.583851612359275, 
              urcrnrlon=1.841589961763497, 
              urcrnrlat=41.598674173123)

ds = gdal.Open("../sample_files/dem.tiff")
elevation = ds.ReadAsArray()

map.imshow(plt.imread('../sample_files/orthophoto.jpg'))

map.imshow(elevation, cmap = plt.get_cmap('terrain'), alpha = 0.5)

plt.show()

plt.read 可以读取任何常规图像并存储到 numpy 数组中
使用 terrian colormap 绘制 elevation （海拔）数据，并且设置 alpha （透明度），从而可以一次看到两个图层
注意： extent 并没有设置，这是因为地图的范围和图像的范围是相同的

第二个例子展示了如何直接在地图上或是在新建的轴上添加 logo：

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
fig = plt.figure()

map = Basemap(projection='ortho', 
              lat_0=0, lon_0=0)

map.drawlsmask(land_color = "#ddaa66", 
               ocean_color="#7777ff",
               resolution = 'l')

x0, y0 = map(1., 31.)
x1, y1 = map(15., 39.)

plt.imshow(plt.imread('../sample_files/by.png'),  extent = (x0, x1, y0, y1))
        
axicon = fig.add_axes([0.1, 0., 0.15, 0.15])
axicon.imshow(plt.imread('../sample_files/by.png'), origin = 'upper')
axicon.set_xticks([])
axicon.set_yticks([])

plt.show()

第一个 logo 使用 extent 参数直接绘制在地图上。坐标被转换为地图单元
第二个 logo 在全球图外部，因此没有地图坐标能对其坐标进行合适的转换

1）轴名为 axicon 的轴被创建，分别确定了轴的 x，y 位置及 width 和 height。这些仅是图形尺寸的一部分。

2）图形被绘制，同时设置 origin 参数，以防止图像上下颠倒

3） xtick 和 ytick 被设置为 null，因此 logo 中并未显示

我从以下例子中学习到了这些。[注8]

pcolor

此函数的行为和 pcolormesh 几乎完全相同。

pcolormesh

创建伪彩色图

pcolormesh(x, y, data, *args, **kwargs)

见 Basemap 系列教程：基本函数绘制栅格数据 部分

plot

在地图上绘制 marker 或 lines。

plot(x, y, *args, **kwargs)

x 和 y 可以是表示投影单元下的marker位置的浮点数，或是绘制 lines 时的点列表
如果 latlon 设置为 True，x 和 y 将被解释为经纬度。旧版本将不会生效。
默认情况下，marker是一个点。[注9]
默认情况下，color是 black（k）。[注10]

第一个例子展示了单个点：

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

map = Basemap(projection='ortho', 
              lat_0=0, lon_0=0)

map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='coral',lake_color='aqua')
map.drawcoastlines()

x, y = map(0, 0)

map.plot(x, y, marker='D',color='m')

plt.show()

第二个例子：如果参数是一个数组，那么输出将是一条线（此例中没有 marker）

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

map = Basemap(projection='ortho', 
              lat_0=0, lon_0=0)

map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='coral',lake_color='aqua')
map.drawcoastlines()

lons = [-10, -20, -25, -10, 0, 10]
lats = [40, 30, 10, 0, 0, -5]

x, y = map(lons, lats)

map.plot(x, y, marker=None,color='m')

plt.show()

quiver

绘制向量场。类似 barbs 方法。

quiver(x, y, u, v, *args, **kwargs)

matplotlib 中的文档介绍更加详细。[注11]

x 和 y 是给定的格点数据位置。如果 latlon 设置为True，这些值将被认为是在地理学坐标系中，否则为地图坐标系。
u 和 v 是左右和上下方向的大小 u 和 v 是以 knot 为单位的左右和上下方向的风的大小。如果投影为非圆柱投影（即除了 cyl，merc，gall，mill [注2]），应该使用 rotate_vector 或 transform_scalar 对 u 和 v 进行旋转

同时还有很多其它可选参数 [注3]。一些常用参数如下：

第五个参数，用于设置箭头的颜色
scale ：控制箭头的长短
minshaft ：值要大于1
pivot 可以改变箭头的旋转点。默认是 ‘tip’，可以改为 ‘middle’

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from osgeo import gdal
import numpy as np


map = Basemap(llcrnrlon=-93.7, llcrnrlat=28., urcrnrlon=-66.1, urcrnrlat=39.5,
              projection='lcc', lat_1=30., lat_2=60., lat_0=34.83158, lon_0=-98.)

ds = gdal.Open("../sample_files/wrf.tiff")
lons = ds.GetRasterBand(4).ReadAsArray()
lats = ds.GetRasterBand(5).ReadAsArray()
u10 = ds.GetRasterBand(1).ReadAsArray()
v10 = ds.GetRasterBand(2).ReadAsArray()
speed = np.sqrt(u10*u10 + v10*v10)

x, y = map(lons, lats)

yy = np.arange(0, y.shape[0], 4)
xx = np.arange(0, x.shape[1], 4)

points = np.meshgrid(yy, xx)

map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='#cc9955', lake_color='aqua', zorder = 0)
map.drawcoastlines(color = '0.15')

map.quiver(x[points], y[points], 
    u10[points], v10[points], speed[points],
    cmap=plt.cm.autumn)

plt.show()

使用 barbs时最主要的问题是点的密度可能会很大，而这一方法又不能自动忽略这些点 1）可仅使用矩阵数据的 1/4 进行绘图 2）矩阵中的点要含有需要的数据点 3）传递给 barbs 的参数可使用点阵进行选取，因此可每 4 个元素选择 1 个值
u 和 v 被认为是地图投影坐标，而不是经纬度坐标。如果是地理学坐标系的话，可使用 rotate_vector 方法进行适当的旋转
计算风速然后设置为 quiver 方法的 color，其数组长度应等于 x，y，u 和 v

scatter

在地图上绘制多个 marker

scatter(x, y, *args, **kwargs)

x 和 y 是要添加到地图上的点的列表
如果 latlon 设置为 True，x 和 y 将被解释为经纬度。旧版本将不会生效。
默认情况下，marker是一个点。[注9]
默认情况下，color是 black（k）。[注10]

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

map = Basemap(projection='ortho', 
              lat_0=0, lon_0=0)

map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='coral',lake_color='aqua')
map.drawcoastlines()

lons = [0, 10, -20, -20]
lats = [0, -10, 40, -20]

x, y = map(lons, lats)

map.scatter(x, y, marker='D',color='m')

plt.show()

streamplot

从向量场绘制流线图。

streamplot(x, y, u, v, *args, **kwargs)

x 和 y 是数组大小和 u， v 数据相同的矩阵，包含在地图坐标中元素的位置 1） matplotlib 文档中解释：x 和 y 必须是均匀间隔的。这就意味着当原始数据源于不同的投影时，data 矩阵必须被重新投影，而且 x 和 y 矩阵也要重新计算，就如此例所示 [注12] 2） makegrid 方法可以计算均匀间隔的网格。最好使用 returnxy = True 属性来获取地图投影单元的网格
u 和 v 是左右和上下方向对应的大小 1） u 和 v 是以 knot 为单位的左右和上下方向的风的大小。如果投影为非圆柱投影（即除了 cyl，merc，gall，mill [注2]），应该使用 rotate_vector 或 transform_scalar 对 u 和 v 进行旋转 2）维度应和 x 和 y 相同
streamplot 可以设置为相同的颜色，或是根据数据设置 1）如果值为标量，所有的 streamplot 都为同一种颜色（取决于colormap） 2）如果数据为和 data 大小相同的数组，那么颜色将随数据变化
cmap 可以设置 colormap
linewidth 用于设置 streamplot 的宽，行为类似 color参数 1）如果为标量，所有的 streamplot 将为固定的宽度 2）如果为数组， streamplot的宽度将随数据变化
density 设置流线的密度
norm 归一化比例，以此设置亮度
arrowsize 控制箭头大小
arrowstyle 设置箭头类型。[注13]
minlength 设置地图坐标中 streamplot 最小长度

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from osgeo import gdal
import numpy


map = Basemap(llcrnrlon=-93.7, llcrnrlat=28., urcrnrlon=-66.1, urcrnrlat=39.5,
              projection='lcc', lat_1=30., lat_2=60., lat_0=34.83158, lon_0=-98.)

ds = gdal.Open("../sample_files/wrf.tiff")
lons = ds.GetRasterBand(4).ReadAsArray()
lats = ds.GetRasterBand(5).ReadAsArray()
u10 = ds.GetRasterBand(1).ReadAsArray()
v10 = ds.GetRasterBand(2).ReadAsArray()
speed = numpy.sqrt(u10*u10 + v10*v10)
xx, yy = map.makegrid(v10.shape[1], v10.shape[0], returnxy=True)[2:4]

map.streamplot(xx, yy, u10, v10, color=speed, cmap=plt.cm.autumn, linewidth=0.5*speed)

map.drawcoastlines(color = '0.15')

plt.show()

注意：

matplotlib 和 basemap 版本要足够高才能使用此函数。如不支持，请升级 matplotlib 和 basemap 版本。

text

添加文本到地图。

text(x, y, s, fontdict=None, withdash=False, **kwargs)

text 方法不属于 Basemap，而是属于 matplotlib。因此，必须使用 plot 或 axis 实例调用
x 和 y 是地图投影坐标。不接受坐标数组，因此，如果要添加多个 label 时，需要多次调用此方法
s 是文本字符串
withdash 设为True 时，将使用 dash‘ 创建字符串
fontdict 用于设置文本性质

text 有很多可选参数：

fontsize
fontweight
ha （horizontal align），比如 center，left 或 right
va （vertical align），比如 center， top 或 bottom
color
bbox 用于在文本周围创建封闭框。比如 bbox=dict(facecolor=’red’, alpha=0.5)

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

map = Basemap(projection='ortho', 
              lat_0=0, lon_0=0)

map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='#cc9955',lake_color='aqua')
map.drawcoastlines()

lon = 3.4
lat = 3.4

x, y = map(lon, lat)

plt.text(x, y, 'Lagos',fontsize=12,fontweight='bold',
                    ha='left',va='bottom',color='k')

lon = 2.1
lat = 41.

x, y = map(lon, lat)

plt.text(x, y, 'Barcelona',fontsize=12,fontweight='bold',
                    ha='left',va='center',color='k',
                    bbox=dict(facecolor='b', alpha=0.2))
plt.show()