文章/答案/技术大牛

发布

Basemap工具函数

共 4 篇文章

Basemap工具函数(2)

Basemap工具函数(1)

Basemap工具函数(4)

Basemap工具函数(3)

清单首页Basemap工具函数文章详情

清单「Basemap工具函数」 03/04

Basemap工具函数(4)

bugsuse

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

注意：

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，这也是最后在地图上所能看到点数
绘制原始数据和插值后的数据