Basemap工具函数(4)

tissot

Tissot 指示图或 Tissot 歪曲椭圆是在地图上显示圆，展示了这些圆是如何适应投影的（即，在不同的位置出现了球面相同的曲率）。通常，不同的位置会出现不同的扭曲度。

tissot(lon_0, lat_0, radius_deg, npts, ax=None, **kwargs)

• lon_0 和 lat_0 是 Tissot 椭圆的位置
• radius_deg 表示多边形的半径
• npts 是多边形的定点数。值越大越接近椭圆

注意：

如果在地图的边缘，圆被分割了（比如从 -179 到 179），此方法不会很好的解决此问题。

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()
for lon in range(0, 360, 20):
    for lat in range(-60, 90, 30):
        map.tissot(lon, lat, 4, 50)
plt.show()

Orthographic 投影的 Tissot 图

Mercator 投影的 Tissot 图

Albers Equal Area 投影的 Tissot 图

transform_scalar

给一个 cylindrical 投影中的标量矩阵及经纬度坐标点，插值这些点到新的矩阵中。

transform_scalar(datin, lons, lats, nx, ny, returnxy=False, checkbounds=False, order=1, masked=False)

• datain 是一个含有标量的二维 numpy 数组
• lons 和 lats 是对应 uin 和 vin 矩阵点的一维 numpy 数组（地理学坐标）。输入的 lon-lat 网格必须是规则的（cyl, merc, mill, cea 和 gall 投影）
• nx 和 ny 是输出网格的ｘ和ｙ的维度。输出网格覆盖了地图，而不是其域外的原点。因此在地图上最终显示的点数是　nx X ny
• returnxy 使此方法返回重新投影后的 lon 和 lat 矩阵。就像调用 basemap 实例一样
• checkbounds 如果为True，将检查 xin 和 yin 是否在 xout 和 yout 边界内。如果为 False，输出数组中那些边界外的值将被裁剪
• masked 如果为True，新网格外的点将被 mask 或置为任意给定值
• order 是插值方法 0 表示最邻近插值；1 表示双线性插值；3 表示三次样条插值，需要scipy.ndimage

注意：

当输入矩阵是不规则的 longitude-latitude 网格（例如，不是 cylindric 投影），此方法将不能正确使用，因为longitude-latitude 网格不是规则的。解决方法见 interp

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
fig=plt.figure(figsize=(9, 3))
map = Basemap(projection='merc', 
              lat_0=0,
              resolution='h',
              llcrnrlon=32, 
              llcrnrlat=31, 
              urcrnrlon=33, 
              urcrnrlat=32)
ds = gdal.Open("../sample_files/dem.tiff")
data = ds.ReadAsArray()
lons = np.linspace(0, 40, data.shape[1])
lats = np.linspace(0, 35, data.shape[0])
ax = fig.add_subplot(121)
ax.set_title('Without transform_scalar')
llons, llats = np.meshgrid(lons, lats)
x, y = map(llons, llats)
map.pcolormesh(x, y, data, vmin=900, vmax=950)
map.drawcoastlines()
ax = fig.add_subplot(122)
ax.set_title('Applying transform_scalar')
data_interp, x, y = map.transform_scalar(data, lons, lats, 40, 40, returnxy=True)
map.pcolormesh(x, y, data_interp, vmin=900, vmax=950)
map.drawcoastlines()
plt.show()

• 本例采用的数据是其它投影和区域的 DEM 数据，因此，我们制作经纬度数据以便使用这些数据
• 使用 linspace 创建等间距的经纬度数组。为了使用 transform_scalar，而且必须是一维数组，因此投影必须是 cylindrical (projections cyl, merc, mill, cea 和 gall)
• 在第一幅图上绘制原始数据

1） 使用 meshgrid 转换经纬度为 二维数组

2） 使用 basemap 实例转换经纬度为 mercator 投影

3） 使用 pcolormesh 绘制结果，并且设置最大值和最小值，因此两幅图的结果是相同的

• 使用 transform_scalar 1） 设置 returnxy 为 True，因此很容易获取新网格的位置 2） 新网格的大小是 40 X 40，因此像素块仍然可见，但是很小。值越大，像素块越小
• 均使用 pcolormesh 绘图，而且最大值和最小值参数设置相同，如果不是的话，命令将根据数据设定值，如果这样的话，颜色看起来可能会不同

注意：

mask似乎并没有效果。

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
fig=plt.figure(figsize=(9, 3))
map = Basemap(projection='robin', 
              lat_0=0, lon_0=0)
lons = np.arange(-180, 190, 60)
lats = np.arange(-90, 100, 30)
data = np.indices((lats.shape[0], lons.shape[0]))
data = data[0] + data[1]
llons, llats = np.meshgrid(lons, lats)
ax = fig.add_subplot(121)
ax.set_title('Without transform_scalar')
x, y = map(llons, llats)
map.pcolormesh(x, y, data, cmap='Paired')
map.drawcoastlines()
ax = fig.add_subplot(122)
ax.set_title('Applying transform_scalar')
data_interp, x, y = map.transform_scalar(data, lons, lats, 50, 30, returnxy=True, masked=True)
map.pcolormesh(x, y, data_interp, cmap='Paired')
map.drawcoastlines()
plt.show()

• 此例使用的数据和 shiftdata 例子中使用的数据相同
• 因为地图覆盖了全球，因此部分输出数组的网格点在地图外 使用 masked = True，这些点将不会有数据，但似乎并没有生效，而且这些点仍然被绘制了，还出现了一些奇怪的现象。

transform_vector

给定向量场的 东西 和 南北 方向分量以及经纬度点，然后对向量进行旋转，使向量场在地图投影上以适当的方向显示。

一些函数（比如 barbs，quiver，streamplot）使用的是向量数据，要求向量分量是地图坐标系（比如 u 是左右方向，v 是上下方向）。如果可用数据是地理学坐标系的（比如东西方向，南北方向），这些坐标必须进行转换，否则所绘制的向量方向会很怪异。这就是 rotate_vertor 方法的目的。

当绘制 barbs，quiver，streamline 时，可用点很少时，需要插值这些点到一个新的矩阵中，从而获取更多的元素以进行绘图。

rotate_vector 方法也能完成同样的工作，但并 没有对点进行插值 。

transform_vector(uin, vin, lons, lats, nx, ny, returnxy=False, checkbounds=False, order=1, masked=False)

• uin 和 vin 是输入数据矩阵。方向是地理学方向，即 u 分量为东西方向， v 分量为 南北方向
• lons 和 lats 是对应 uin 和 vin 矩阵点的一维 numpy 数组（地理学坐标）。输入的 lon-lat 网格必须是规则的（cyl, merc, mill, cea 和 gall 投影）
• nx 和 ny 是输出网格的ｘ和ｙ的维度。输出网格覆盖了地图，而不是其域外的原点。因此在地图上最终显示的点数是　nx X ny
• returnxy 使此方法返回重新投影后的 lon 和 lat 矩阵。就像调用 basemap 实例一样
• checkbounds 如果为True，将检查 xin 和 yin 是否在 xout 和 yout 边界内。如果为 False，输出数组中那些边界外的值将被裁剪
• masked 如果为True，新网格外的点将被 mask 或置为任意给定值
• order 是插值方法 0 表示最邻近插值；1 表示双线性插值；3 表示三次样条插值，需要scipy.ndimage

  当输入矩阵是不规则的 longitude-latitude 网格（例如，不是 cylindric 投影），此方法将不能正确使用，因为longitude-latitude 网格不是规则的。解决方法见 interp

  from mpl_toolkits.basemap import Basemap
  import matplotlib.pyplot as plt
  from osgeo import gdal
  import numpy as np
  map = Basemap(projection='sinu', 
                lat_0=0, lon_0=0)
  lons = np.linspace(-180, 180, 8)
  lats = np.linspace(-90, 90, 8)
  v10 = np.ones((lons.shape)) * 15
  u10 = np.zeros((lons.shape))
  u10, v10 = np.meshgrid(u10, v10)
  u10_rot, v10_rot, x_rot, y_rot = map.transform_vector(u10, v10, 
                                                        lons, lats, 
                                                        15, 15, 
                                                        returnxy=True)
  map.drawmapboundary(fill_color='#9999FF')
  map.fillcontinents(color='#ddaa66', lake_color='#9999FF', zorder = 0)
  map.drawcoastlines(color = '0.15')
  #Drawing the original points
  lons, lats = np.meshgrid(lons, lats)
  x, y = map(lons, lats)
  map.barbs(x, y, u10, v10, 
      pivot='middle', barbcolor='#333333')
  #Drawing the rotated & interpolated points
  map.barbs(x_rot, y_rot, u10_rot, v10_rot, 
      pivot='middle', barbcolor='#ff5555')
  plt.show()

  复制
  • lons 和 lats 是覆盖全球的等间距网格
  • v10 和 u10 创建后，呈现为南北风向（v10 = 10， u10 = 0）.使用 meshgrid 转换为 二维数组
  • 一旦数据被创建了，可以使用 transform_vector 旋转和插值向量并返回新的网格 设置 nx 和 ny 为15，在地图投影上新的网格将是 15 x 15，这也是最后在地图上所能看到点数
  • 绘制原始数据和插值后的数据