使用plot绘图

在地图上绘制标记或线条。

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

第一个示例显示单个点:

In [1]:
%matplotlib inline
import warnings
warnings.filterwarnings('ignore')
In [2]:
import os
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt 
import numpy as np
my_map = Basemap(projection='merc', lat_0=57, lon_0=-135,
                 resolution='h', area_thresh=0.1,
                 llcrnrlon=-136.25, llcrnrlat=56.0,
                 urcrnrlon=-134.25, urcrnrlat=57.75)
In [3]:
my_map.drawcoastlines()
my_map.drawcountries()
my_map.fillcontinents(color='coral')
my_map.drawmapboundary()
lon = -135.3318
lat = 57.0799
x, y = my_map(lon, lat)
my_map.plot(x, y, 'bo', markersize=12)

plt.show()

x和y可以是带有投影单元中标记位置的浮点,也可以是带有绘制线的点的列表。

如果lat与lon关键字设置为True,则x,y将解释为经度和纬度,以度为单位。不会在旧的Basemap版本中工作。

  • 默认情况下,标记是一个点。此页面解释所有选项。
  • 默认情况下,颜色为black(k)。

现在让我们给地图添加一些点。例如在Sitan,Baranof岛上最大的社区,我们添加一个点,显示Sitka的位置。在plt.show()之前添加以下行:

In [4]:
my_map = Basemap(projection='merc', lat_0=57, lon_0=-135,
                 resolution='h', area_thresh=0.1,
                 llcrnrlon=-136.25, llcrnrlat=56.0,
                 urcrnrlon=-134.25, urcrnrlat=57.75)

my_map.drawcoastlines()
my_map.drawcountries()
my_map.fillcontinents(color='coral')
my_map.drawmapboundary()

lons = [-135.3318, -134.8331, -134.6572]
lats = [57.0799, 57.0894, 56.2399]
x, y = my_map(lons, lats)
my_map.plot(x, y, 'bo', markersize=10)

plt.show()

在地图中绘制点中,通常使用plot方法:

  • 当您具有点的经度和纬度时,使用Basemap实例计算地图坐标中点的位置
  • 如果lat和lon的关键字设置为True,则x,y将解释为经度和纬度,以度为单位。不会在旧Basemap的版本中工作.

绘图需要地图坐标中的x和y位置,在标记和颜色都默认情况下,标记是一个点。 此页面解释所有选项。默认情况下,颜色为black(k)。

这里唯一不明显的是bo参数,它告诉底图对点使用蓝色圆圈。有相当多的颜色和符号,你可以使用。默认标记大小为6,但在此特定地图上太小。12的markersize显示在这张地图上。

绘制单个点的方法是很好的,但我们经常想在地图上绘制一个大的点集。以Baranof岛上的两个其他社区为例,展示这两个社区在这个地图上的位置。我们将点的纬度和经度存储在两个单独的列表中,将它们映射到x和y坐标,并在地图上绘制这些点。因为会在地图上还有更多的点,我们还想应该要稍微减小标记大小:

标签点

现在让我们标记这三点。首先创建一个标签列表,并循环该列表。我们需要在此循环中为包括x和y值的每个点,因此Basemap可以确定放置每个标签的位置。

In [5]:
my_map = Basemap(projection='merc', lat_0=57, lon_0=-135,
                 resolution='h', area_thresh=0.1,
                 llcrnrlon=-136.25, llcrnrlat=56.0,
                 urcrnrlon=-134.25, urcrnrlat=57.75)

my_map.drawcoastlines()
my_map.drawcountries()
my_map.fillcontinents(color='coral')
my_map.drawmapboundary()

lons = [-135.3318, -134.8331, -134.6572]
lats = [57.0799, 57.0894, 56.2399]
x, y = my_map(lons, lats)
my_map.plot(x, y, 'bo', markersize=10)
# 标记了三个点。

labels = ['Sitka', 'Baranof Warm Springs', 'Port Alexander']
# 每个点分别赋值。

for label, xpt, ypt in zip(labels, x, y):
    plt.text(xpt, ypt, label)
plt.show()

如果参数是数组,则输出是一行(在这种情况下不带标记):

In [6]:
import numpy as np

my_map = Basemap(projection='merc', lat_0 = 57, lon_0 = -135,
    resolution = 'h', area_thresh = 0.1,
    llcrnrlon=-136.25, llcrnrlat=56.0,
    urcrnrlon=-134.25, urcrnrlat=57.75)

my_map.drawcoastlines()
my_map.drawcountries()
my_map.fillcontinents(color = 'coral')
my_map.drawmapboundary()

lons = [-135.3318, -134.8331, -134.6572]
lats = [57.0799, 57.0894, 56.2399]
x,y = my_map(lons, lats)
my_map.plot(x, y, 'bo', markersize=10)

labels = ['Sitka', 'Baranof Warm Springs', 'Port Alexander']
for label, xpt, ypt in zip(labels, x, y):
    plt.text(xpt+10000, ypt+5000, label)

plt.show()

我们通过上图可以看到的城镇已被标记,但标签开始在点的顶部。而且它们不是正确的点,我们可以向这些点添加偏移量。让我们将所有标签向上和向右移动一点。 (如果你好奇,这些偏移在地图投影坐标中,以米为单位,这意味着我们的代码实际上将标签放置在东部10公里,实际城镇北部5公里。)

In [7]:
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')
# 顺序连接每个点 颜色为m
plt.show()

quiver

quiver(x, y, u, v, *args, **kwargs)
In [8]:
my_map = Basemap(projection='merc', lat_0=57, lon_0=-135,
                 resolution='h', area_thresh=0.1,
                 llcrnrlon=-136.25, llcrnrlat=56.0,
                 urcrnrlon=-134.25, urcrnrlat=57.75)

my_map.drawcoastlines()
my_map.drawcountries()
my_map.fillcontinents(color='coral')
my_map.drawmapboundary()

lons = [-135.3318, -134.8331, -134.6572]
lats = [57.0799, 57.0894, 56.2399]
x, y = my_map(lons, lats)
my_map.plot(x, y, 'bo', markersize=10)

labels = ['Sitka', 'Baranof\n  Warm Springs', 'Port Alexander']
x_offsets = [10000, -20000, -25000]
y_offsets = [5000, -50000, -35000]

for label, xpt, ypt, x_offset, y_offset in zip(labels, x, y, x_offsets, y_offsets):
    plt.text(xpt + x_offset, ypt + y_offset, label)


plt.show()
  • x和y给出了网格数据的位置,如果lat和lon自变量为真,则这些值应该在地理坐标中。如果没有,这些值应该在地图坐标中。
  • u和v是表示左右和上下幅度。

注意,它们不在南北或东西。如果输入投影具有非圆柱形投影(除cyl,merc,cyl,gall 和mill之外的投影),则应使用rotate_vector或transform_scalar方法旋转u和v。

  • 第五个参数(可选)设置一个值将颜色分配给箭头。
  • scale使箭头更长或更短。
  • 低于1的值的箭头增长。
  • 枢轴改变箭头的旋转点。默认情况下是“tip”,但可以更改为“middle”。

这是更好的,但在这个尺度的地图上,相同的偏移量并不适用于所有点。我们可以单独绘制每个标签,但是最好是做出两个列表来存储偏移: 有没有容易的方法来保持“Baranof温泉”越过边界,但在标签中使用换行符使它更清晰一些。现在我们知道如何添加点到地图,我们可以移动到更大的数据集。

多点Scatter

如果图像中有多个点,您可能更喜欢分散方法。首先创建一条连接它们的线将点数组传递到plot方法:

在地图上绘制多个标记。

scatter(x, y, *args, **kwargs)
In [9]:
my_map = Basemap(projection='merc', lat_0=57, lon_0=-135,
                 resolution='h', area_thresh=0.1,
                 llcrnrlon=-136.25, llcrnrlat=56.0,
                 urcrnrlon=-134.25, urcrnrlat=57.75)

my_map.drawcoastlines()
my_map.drawcountries()
my_map.fillcontinents(color='coral')
my_map.drawmapboundary()

lons = [-135.3318, -134.8331, -134.6572]
lats = [57.0799, 57.0894, 56.2399]
x, y = my_map(lons, lats)
my_map.plot(x, y, 'bo', markersize=10)

# labels = ['Sitka', 'Baranof\n  Warm Springs', 'Port Alexander']
x_offsets = [10000, -20000, -25000]
y_offsets = [5000, -50000, -35000]

for label, xpt, ypt, x_offset, y_offset in zip(labels, x, y, x_offsets, y_offsets):
    plt.text(xpt + x_offset, ypt + y_offset, label)

plt.show()
  • x和y是要作为标记添加到地图的点列表。
  • 如果latlon关键字设置为True,则x,y将解释为经度和纬度,以度为单位。不会在旧Basemap版本中工作。
  • 默认情况下,标记是一个点。
  • 默认情况下,颜色为黑色(k)。
In [10]:
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()

请记住,调用Basemap实例可以使用列表来完成,所以坐标变换也会立即完成。

散点法中的格式选项与绘图中的相同。

Streamplot

从矢量字段绘制流线图。

streamplot(x, y, u, v, *args, **kwargs)
In [11]:
my_map = Basemap(projection='merc', lat_0=57, lon_0=-135,
                 resolution='h', area_thresh=0.1,
                 llcrnrlon=-136.25, llcrnrlat=56.0,
                 urcrnrlon=-134.25, urcrnrlat=57.75)

my_map.drawcoastlines()
my_map.drawcountries()
my_map.fillcontinents(color='coral')
my_map.drawmapboundary()

lons = [-135.3318, -134.8331, -134.6572]
lats = [57.0799, 57.0894, 56.2399]
x, y = my_map(lons, lats)
my_map.plot(x, y, 'bo', markersize=10)

labels = ['Sitka', 'Baranof\n  Warm Springs', 'Port Alexander']
x_offsets = [10000, -20000, -25000]
y_offsets = [5000, -50000, -35000]

for label, xpt, ypt, x_offset, y_offset in zip(labels, x, y, x_offsets, y_offsets):
    plt.text(xpt + x_offset, ypt + y_offset, label)


plt.show()
  • x和y是与u和v数据大小相同的矩阵,包含地图坐标中元素的位置。
    • 正如文档解释,x和y必须均匀分布。这意味着当原始值来自不同的投影时,必须重新投影 数据矩阵,并重新计算x和y矩阵,如在示例中所见。
    • 要均匀计算间隔的网格,可以使用方法makegrid。最好使用returnxy = True属性来获取地图投影单元中的网格。
  • u和v是左右和上下幅度

    • 注意,它们不在南北和东西。如果输入的投影具有非圆柱形投影(除cyl,merc,cyl,gall 和mill之外的投影),则应使用rotate_vector或transform_scalar方法旋转u和v

    • 尺寸必须与x和y相同

  • 颜色可以为流线设置相同的颜色,或根据数据进行更改:

    • 如果值是标量,则所有流线将具有指示的颜色,这取决于色彩映射。
    • 如果值是与数据(在示例中的风的模块)相同大小的数组,则颜色将根据它而改变,使用颜色图。
  • cmap设置颜色图
  • lineewidth以类似于颜色的方式设置流线的宽度
    • 如果是标量,则所有流线具有指示的宽度。
    • 如果是一个数组,流线宽度将根据数组的值而改变。
  • 可以通过密度设置流线图的接近程度。 1值意味着域被划分为30×30,每个扇区只有一条流线。如果传递具有两个元素的列表,则将x和y密度设置为不同的值
  • norm标准化标尺以设置亮度数据。
  • arrowize缩放箭头大小。
  • arrowstyle stes箭头样式。文档在FancyArrowPatch minlength中。
  • 设置地图坐标中流线的最小长度。

text

在地图上绘制text。

text(x, y, s, fontdict=None, withdash=False, **kwargs)
  • text方法不属于Basemap,但直接属于matplotlib,因此必须从图或轴实例中调用。
In [12]:
map = Basemap(projection='ortho', lat_0=0, lon_0=105)

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()
  • x和y是地图投影中的坐标。不接受坐标数组,因此要添加多个标签,应多次调用该方法。
  • s是文本字符串。
  • 当设置为true时,withdrawash将创建一个带有短划线的text。
  • fontdict可用于对text属性进行分组。
  • text可以有许多选项,如:
    • fontsize字体大小。
    • fontweight的字体重量,如粗体。
    • ha horitzontal对齐,像中心,左或右。
    • va垂直对齐,如中心,顶部或底部。
    • 图像的颜色。
    • bbox在文本周围创建一个框:bbox = dict(facecolor ='red',alpha = 0.5)。

使用shapefile

Basemap处理矢量文件的方式与其他库非常不同,值得注意。

基本用法

让我们从绘制shapefile最简单的方法开始:

In [13]:
%matplotlib inline
import warnings
warnings.filterwarnings('ignore')
import os
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
  • 第一个参数shapefile名称没有shp扩展名。该库假定所有shp、sbf和shx文件将以此给定名称。
  • 第二个参数是稍后从Basemap实例访问shapefile信息的名称,我们将在后面介绍。
In [16]:
map = Basemap(llcrnrlon=-0.5,llcrnrlat=39.8,urcrnrlon=4.,urcrnrlat=43.,
             resolution='i', projection='tmerc', lat_0 = 39.5, lon_0 = 1)

map = Basemap(projection='robin', lat_0=0, lon_0=-100,
                 resolution='l', area_thresh=1000.0)

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

map.readshapefile('/gdata/GSHHS_c', 'comarques')
plt.show()

有一些限制:

  • 文件必须为EPSG:4326或lat / lon坐标。如果您的文件不是,您可以使用ogr2ogr来转换它
  • 元素只能有2个维度
  • 此示例将仅在元素是多边形或折线时显示
  • 如图所示,结果将是多边形(或折线)的边界。要填充它们,请查看最后一部分填充多边形

读取点数据

绘制点有点复杂。首先,读取shapefile,然后可以使用散点图、绘图点或matplotlib函数来绘制点,以更好地满足需要。

例子显示在暴风雨期间在加泰罗尼亚落下的闪电

第二个参数已命名为lightnings和Basemap实例映射,所以shapefile元素可以使用map.lightnings几何和map.lightnings_info参数访问元素字段

shapefile方法返回具有元素数量的序列,具有此处定义的代码的几何类型和边界框

第17行显示了迭代所有元素的可能方法:

  • zip将使用其字段值连接每个几何
  • 每个元素as当迭代一个dict可以迭代一个for。
  • 在该示例中,称为振幅的场可以用于猜测闪电是具有正电流还是负电流,并且针对每种情况绘制不同的符号
  • 使用方法图绘制点,使用marker属性更改符号

多边形信息

此示例显示如何使用shapefile属性选择一些几何。

填充多边形

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

map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='#ddaa66',lake_color='aqua')
# 用#ddaa66颜色填充大陆
map.drawcoastlines()

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

plt.show()

# plt.savefig( get_tmp_file(__file__, 2))

matplotlib正在使的是用一个称为PatchCollection的类,它是一个集合形状以添加填充多边形,如官方文档中所解释。

在这种情况下,形状是Polygon类型。要创建它,坐标必须是numpy数组。第二个参数设置要关闭的多边形。

绘制栅格数据

想要绘制轮廓线或填充轮廓线(等值线)和pcolor / pcolormesh的栅格,有两种主要方法,一是创建伪彩色图。

In [ ]:
%matplotlib inline
import warnings
warnings.filterwarnings('ignore')
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from osgeo import gdal
from numpy import linspace
from numpy import meshgrid

para = {
    'projection': 'tmerc',
    'lat_0': 0,
    'lon_0': 3,
    'llcrnrlon': 1.819757266426611,
    'llcrnrlat': 41.583851612359275,
    'urcrnrlon': 1.841589961763497,
    'urcrnrlat': 41.598674173123
}
dem_tif = '/gdata/sample_files/dem.tiff'
  • 使用与栅格文件相同的扩展名创建映射,这个方法比较简单。
In [ ]:
p1 = plt.subplot(121)
map = Basemap(**para)

ds = gdal.Open(dem_tif)
data = ds.ReadAsArray()

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

xx, yy = meshgrid(x, y)

cs = map.contour(xx, yy, data, range(400, 1500, 100), cmap=plt.cm.cubehelix)

p2 = plt.subplot(122)

map = Basemap(**para)

ds = gdal.Open(dem_tif)
data = ds.ReadAsArray()

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

xx, yy = meshgrid(x, y)

cs = map.contour(xx, yy, data, range(400, 1500, 100), cmap=plt.cm.cubehelix)
plt.clabel(cs, inline=True, fmt='%1.0f', fontsize=12, colors='k')

# fig_index += 1
# plt.savefig(get_tmp_file(__file__, fig_index))
plt.show()
# plt.clf()
  • 在绘制轮廓之前,必须创建两个矩阵,包含数据矩阵中每个点的x和y坐标的位置。
  • linspace是一个numpy函数,用于创建一个包含n个元素的从初始值到结束值的数组。在这种情况下,地图坐标从0到map.urcrnrx或map.urcrnry,并且具有与数据数组 data.shape [1]和data.shape [0]相同的大小。
  • meshgrid是一个numpy函数,它接受两个数组并与它们创建一个矩阵。这一点很重要,因为x坐标在每一列重复,y就在每一行。
  • contourf方法将获取x,y和数据矩阵,并将它们绘制在默认色图(称为jet)中,并且自动设置级数。
  • 级数可以在数据数组后面定义,您可以在轮廓截面中看到。这可以通过两种方式完成。
  • 指示级数的整数。数据数组的极值将指示色标的极值。
  • 具有每个级别的值的列表。范围函数对于设置它们是有用的,即范围(0,3000,100)以每100个单位设置级别。
  • 数据必须按照contourf情况进行准备:
    • 水平设置使用范围函数。我们正在处理高度,所以从400米到1400米,每100米创建一条等高线
    • 色彩映射不是默认的jet。这是通过cmap参数传递到cubehelix色彩映射来完成的。
    • 标签可以设置为轮廓法(但不能设置为contourf)
    • inline使线下的轮廓线被删除
    • fmt格式化数字
    • fontsize设置标签字体的大小
    • 标签颜色的设置。默认情况下,与轮廓线相同

数据必须按照contourf情况进行准备

可以像在轮廓示例中那样更改色彩映射.

注:pcolor和pcolormesh非常相似。在这里你可以看到一个很好的解释

contour

创建等高线图。

contour(x, y, data)
In [ ]:
# from numpy import meshgrid

p1 = plt.subplot(121)

map = Basemap(**para)

ds = gdal.Open(dem_tif)
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)

创建等高线图。

In [ ]:
# plt.show()


p2 = plt.subplot(122)

map = Basemap(**para)

ds = gdal.Open(dem_tif)
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.pcolormesh(xx, yy, data)

# plt.savefig(get_tmp_file(__file__, fig_index))

plt.show()
# plt.clf()
# plt.close()
  • x和y是与数据相同大小的矩阵,包含地图坐标中元素的位置。
  • 数据是包含要绘制的数据值的矩阵。
  • 可以传递第四个参数,其中包含创建contour时要使用的级别列表。
  • 默认色彩映射是jet,但是参数cmap可用于更改行为。
  • 参数tri = True使得网格被假定为非结构化的。
  • 其他可能的参数记录在matplotlib函数的docs中。
  • 可以在contour结果中添加标签,如在基本功能部分的contour显示。

Contourf

创建填充轮廓图。

contourf(x, y, data)
  • x和y是与数据相同大小的矩阵,包含地图坐标中元素的位置。
  • 数据是包含要绘制的数据值的矩阵。
  • 默认色彩映射是jet,但是参数cmap可用于更改行为。
  • 第四个参数可以传递,包含创建contourf时要使用的级别列表。
  • 参数tri = True使得网格被假定为非结构化的。
  • 其他可能的参数记录在matplotlib函数的docs中。

pcolor

此函数的运行几乎与pcolormesh中的一致。

pcolormesh

创建伪色图。

pcolormesh(x, y, data, *args, **kwargs)
  • x和y是与数据相同大小的矩阵,包含地图坐标中元素的位置。
  • 数据是包含要绘制的数据值的矩阵。
  • 默认色彩映射是jet,但是参数cmap可用于更改行为。
  • 其他可能的参数记录在matplotlib函数的docs中。

使用shapefile剪切栅格

1、获取一些数据

该示例绘制取自SRTM的一些高程数据。在寻找一些选项后,最容易使用的是这一个:http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp 或直接下载文件:http://srtm.csi.cgiar.org/SRT-ZIP/SRTM_V41/SRTM_Data_GeoTiff/srtm_37_04.zip

shapefile将是安道尔的边界,取自自然地球

这是一个分辨率低,但运行良好的示例。

2、代码

示例脚本使用pyshp读取shapefile。当然,ogr也可以使用,但不能是fiona,因为失败时与gdal在同一个脚本。

In [ ]:
from matplotlib.path import Path
from matplotlib.patches import PathPatch
from matplotlib import ReadAsArray
import  shapefile
import numpy
from geopandas import GeoDataFrame

fig = plt.figure()
ax = fig.add_subplot(111)

sf = shapefile.Reader("../gdata/ne_110m_admin_0_countries.shp")

for shape_rec in sf.shapeRecords():
    if shape_rec.record[3] == 'Andorra':
        vertices = []
        codes = []
        pts = shape_rec.shape.points
        prt = list(shape_rec.shape.parts) + [len(pts)]
        for i in range(len(prt) - 1):
            for j in range(prt[i], prt[i+1]):
                vertices.append((pts[j][0], pts[j][1]))
            codes += [Path.MOVETO]
            codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)
            codes += [Path.CLOSEPOLY]
        clip = Path(vertices, codes)
        clip = PathPatch(clip, transform=ax.transData)


m = Basemap(llcrnrlon=1.4,
    llcrnrlat=42.4,
    urcrnrlon=1.77,
    urcrnrlat=42.7,
    resolution = None,
    projection = 'cyl')

ds = gdal.Open('../srtm_37_04.tif')
data = ds.ReadAsArray()

gt = ds.GetGeoTransform()
x = numpy.linspace(gt[0], gt[0] + gt[1] * data.shape[1], data.shape[1])
y = numpy.linspace(gt[3], gt[3] + gt[5] * data.shape[0], data.shape[0])

xx, yy = numpy.meshgrid(x, y)

cs = m.contourf(xx,yy,data,range(0, 3600, 200))

for contour in cs.collections:
        contour.set_clip_path(clip)

plt.show()

①、读取shapefile和剪切

要剪切图像,需要一个底图路径。代码中突出显示的行是执行该作业的行。

Matplotlib路径由两个数组组成。一个具有点(在脚本中称为顶点),另一个具有每个点的函数(称为代码)。

在我们的例子中,只需要使用直线,所以将有一个MOVETO来指示多边形的开始,许多LINETO创建段和一个CLOSEPOLY关闭它。

当然,只有安道尔的多边形必须使用。我可以从shapefile属性中得到它。

prt数组用于管理多边形集合,情况并非如此,但是代码将为多边形集合创建正确的裁剪。

该路径使用Path函数创建,然后添加到PathPatch,以便能够将其用作封闭多边形。注意trasnform = ax.transData属性。这假设多边形坐标是数据中使用的那些(在我们的例子中是经度和纬度)。

剪切本身在行48和49中进行。对于每个绘制元素,都使用set_clip_path方法,其擦除剪切对象外部的所有部分。

②、绘制地图

地图是照常绘制的。我使用了一个latlon投影,所以栅格和shapefile的所有值都可以直接使用。如果输出栅格在另一个投影中,则应使用输出投影(m(pts [j] [0],pts [j] [1]))将shapefile坐标附加到路径。

x和y坐标从GDAL地理变换计算,然后变成使用meshgrid的矩阵。

绘制标注

annotate

创建带有指示关注点的箭头的文本。

annotate(args, *kwargs)    
In [ ]:
%matplotlib inline
import warnings
warnings.filterwarnings('ignore')
In [ ]:
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
  • 文本方法不属于Basemap,但属于matplotlib,因此必须从图或轴实例中调用。
In [ ]:
fig = plt.figure(figsize=(6, 4)) 
plt.subplots_adjust(left=0.05, right=0.95, top=0.90, bottom=0.05, wspace=0.15, hspace=0.05)
ax = plt.subplot(111)

m = Basemap(resolution='i', projection='merc', llcrnrlat=10.0, urcrnrlat=55.0, llcrnrlon=60., urcrnrlon=140.0)
m.drawcountries(linewidth=0.5)
m.drawcoastlines(linewidth=0.5)

m.drawparallels(np.arange(10., 55., 10.), labels=[1, 0, 0, 0], color='black', labelstyle='+/-', linewidth=0.2,
                dashes=(None, None))  # draw parallels,dashes=[1,0],
m.drawmeridians(np.arange(60., 140., 10.), labels=[0, 0, 0, 1], color='black', labelstyle='+/-', linewidth=0.2,
                dashes=(None, None))  # draw meridians ,dashes=[1,0]
x, y = m(116.4204, 40.21244)  # Bejing
x2, y2 = m(125.27538, 43.83453)  # Changchun
plt.annotate('Beijing', xy=(x, y), xycoords='data',
             # xytext=(x2, y2), textcoords='offset points',
             xytext=(x2, y2), textcoords='data',
             color='r',
             arrowprops=dict(arrowstyle="fancy", color='g')
             )
plt.show()

上面第一行的第一个参数是文本字符串 。用于对自定义子图的调整。

  • xy是具有由箭头指向的点的x和y坐标的列表。这将是interprete依赖于xycoords参数 。xycoords表示xy中使用的坐标类型:
  • 数据意味着数据使用的坐标(投影坐标)。
  • 偏移点表示点的偏移 。
  • axis pixels表示从轴的左下角开始的像素。
  • 其他选项位于注释文档。
  • text是与箭头一起的文字,xy是箭头所在的位置。
  • textcoords表示xytext中使用的坐标类型,具有与xycoords中相同的选项 。
  • arrowprops此可指定箭头的属性,如Line2D文档中所述 。
  • 颜色文本的颜色。