目录

上一个主题

9.3. 使用plot绘图

下一个主题

9.5. 球面距离案例

>>> from helper import info; info()
页面更新时间: 2021-08-04 20:33:14
运行环境:
     Linux发行版本: Debian GNU/Linux 11 (bullseye)
     操作系统内核: Linux-5.10.0-8-amd64-x86_64-with-glibc2.31
     Python版本: 3.9.2

9.4. 使用数据

9.4.1. 使用shapefile

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

基本用法

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

>>> %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信息的名称,我们将在后面介绍。

>>> 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()
>>>
_images/basemap-data_5_0.png

有一些限制:

  • 文件必须为EPSG:4326或lat / lon坐标。如果您的文件不是,您可以使用ogr2ogr来转换它

  • 元素只能有2个维度

  • 此示例将仅在元素是多边形或折线时显示

  • 如图所示,结果将是多边形(或折线)的边界。要填充它们,请查看最后一部分填充多边形

9.4.2. 读取点数据

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

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

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

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

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

  • zip将使用其字段值连接每个几何

  • 每个元素as当迭代一个dict可以迭代一个for。

  • 在该示例中,称为振幅的场可以用于猜测闪电是具有正电流还是负电流,并且针对每种情况绘制不同的符号

  • 使用方法图绘制点,使用marker属性更改符号

9.4.3. 多边形信息

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

  • 要迭代所有元素,请像上一个示例一样使用zip

  • 有一个名为’nombre’的字段可以用来过滤。在该示例中,仅选择值“Selva”即可。

  • 要绘制线,x和y坐标必须是分离的数组,但几何将每个点作为一对。有一个解释如何做到网址http://stackoverflow.com/questions/13635032/what-is-the-inverse-function-of-zip-in-python’_

  • 使用plot方法绘制形状,消除标记以获得一条线

9.4.4. 绘制栅格数据

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

>>> %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'
  • 使用与栅格文件相同的扩展名创建映射,这个方法比较简单。

>>> 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')
>>>
>>>
>>>
>>> plt.show()
---------------------------------------------------------------------------

ProjError                                 Traceback (most recent call last)

<ipython-input-5-d0bedafc6048> in <module>
      1 p1 = plt.subplot(121)
----> 2 map = Basemap(**para)
      3
      4 ds = gdal.Open(dem_tif)
      5 data = ds.ReadAsArray()


/usr/lib/python3/dist-packages/mpl_toolkits/basemap/__init__.py in __init__(self, llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat, llcrnrx, llcrnry, urcrnrx, urcrnry, width, height, projection, resolution, area_thresh, rsphere, ellps, lat_ts, lat_1, lat_2, lat_0, lon_0, lon_1, lon_2, o_lon_p, o_lat_p, k_0, no_rot, suppress_ticks, satellite_height, boundinglat, fix_aspect, anchor, celestial, round, epsg, ax)
   1129             # replace coastsegs with line segments (instead of polygons)
   1130             self.coastsegs, types =\
-> 1131             self._readboundarydata('gshhs',as_polygons=False)
   1132         # create geos Polygon structures for land areas.
   1133         # currently only used in is_land method.


/usr/lib/python3/dist-packages/mpl_toolkits/basemap/__init__.py in _readboundarydata(self, name, as_polygons)
   1410                                     # transformation from lat/lon to
   1411                                     # map projection coordinates.
-> 1412                                     bx, by = self(blons, blats)
   1413                                     if not as_polygons or len(bx) > 4:
   1414                                         polygons.append(list(zip(bx,by)))


/usr/lib/python3/dist-packages/mpl_toolkits/basemap/__init__.py in __call__(self, x, y, inverse)
   1190             except TypeError:
   1191                 y = [_dg2rad*yy for yy in y]
-> 1192         xout,yout = self.projtran(x,y,inverse=inverse)
   1193         if self.celestial and inverse:
   1194             try:


/usr/lib/python3/dist-packages/mpl_toolkits/basemap/proj.py in __call__(self, *args, **kw)
    291             outxy = self._proj4(xy, inverse=inverse)
    292         else:
--> 293             outx,outy = self._proj4(x, y, inverse=inverse)
    294         if inverse:
    295             if self.projection in ['merc','mill','gall']:


/usr/lib/python3/dist-packages/pyproj/proj.py in __call__(self, longitude, latitude, inverse, errcheck, radians)
    180         else:
    181             direction = TransformDirection.FORWARD
--> 182         return self.transform(
    183             xx=longitude,
    184             yy=latitude,


/usr/lib/python3/dist-packages/pyproj/transformer.py in transform(self, xx, yy, zz, tt, radians, errcheck, direction)
    486             intime = None
    487         # call pj_transform.  inx,iny,inz buffers modified in place.
--> 488         self._transformer._transform(
    489             inx,
    490             iny,


pyproj/_transformer.pyx in pyproj._transformer._Transformer._transform()


ProjError: x, y, z, and time must be same size
_images/basemap-data_11_1.png
  • 在绘制轮廓之前,必须创建两个矩阵,包含数据矩阵中每个点的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)
>>>
>>> # 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)
>>>

创建等高线图。

>>> # 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()
  • 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中。

9.4.5. 绘制标注

annotate

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

annotate(args, *kwargs)
>>> %matplotlib inline
>>> import warnings
>>> warnings.filterwarnings('ignore')
>>> from mpl_toolkits.basemap import Basemap
>>> import matplotlib.pyplot as plt
>>> import numpy as np

文本方法不属于Basemap,但属于 Matplotlib,因此必须从图或轴实例中调用。

>>> 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文档中所述 。

  • 颜色文本的颜色。