Basemap工具函数(2)

bugsuse

发布于 2020-04-21 17:31:36

1.6K0

发布于 2020-04-21 17:31:36

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

interp

格点插值。

对于创建平滑图形或使用 barbs 或 quiver 绘图时非常有用。当使用 maskoceans 函数时也非常有用。

interp(datain, xin, yin, xout, yout, checkbounds=False, masked=False, order=1)

不是 Basemap 实例方法，而是独立函数
datain 是要插值的输入数组。可以是二维数组
xin 和 yin 是要插值数据的坐标数组，均为一维数组 1）使用的输入网格点必须是规则的 2） y 必须是递增的
xout 和 yout 是输出数组的坐标数组。必须是二维的 numpy 数组
checkbounds 如果为 True，将会检查 xin 和 yin 中的值是否在 xout 和 yout 边界内，如果为False，对于那些边界外的值，输出数组的值将被剪切为边界的值
masked 如果为 True，新网格外的点将被 mask，或置为给定的任意值
order 设置插值方法

0 表示最邻近插值，1 表示双线性插值， 2 表示三次样条插值，需要安装 scipy.ndimage

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

map = Basemap(llcrnrlon=-82., llcrnrlat=28., urcrnrlon=-79., urcrnrlat=29.5,
              projection='lcc', lat_1=30., lat_2=60., lat_0=34.83158, lon_0=-98.,
              resolution='i')

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)

x2 = np.linspace(x[0][0],x[0][-1],x.shape[1]*2)
y2 = np.linspace(y[0][0],y[-1][0],y.shape[0]*2)

x2, y2 = np.meshgrid(x2, y2)

u10_2 = interp(u10,  x[0], np.flipud(y[:, 0]), x2, np.flipud(y2),order=1)
v10_2 = interp(v10,  x[0], np.flipud(y[:, 0]), x2, np.flipud(y2),order=1)

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

map.barbs(x, y, u10, v10, 
    pivot='middle', barbcolor='#555555')

map.barbs(x2, y2, u10_2, v10_2, 
    pivot='middle', barbcolor='#FF3333')

plt.show()

译注：以下是关于上述代码的部分解释

采用 barb 中的例子
使用 basemap 实例计算 wind barb(x,y) 的位置
使用 linspace 创建新的网格 1）因为网格在地图投影中，x 和 y 是二维数组，需要转换为一维数组

2）新的网格被创建，其元素数是之前的2倍，但拥有相同的边界

3）新的网格必须是二维数组。上例中使用 meshgrid 方法创建

调用 interp 1）因为 x 和 y 是二维数组，因此只能传递一列给 interp 2） y 轴不是单调递增的（坐标是由北纬到南纬），因此必须使用 numpy.flipud 方法反转。输出必须要在此反转才能有合适的顺序
一旦新的网格被创建，就可以使用 barbs 绘图了

灰色为原始点，红色为插值点。

is_land

如果所选点在 land 则返回 True，如果在 ocean 或 lake 则返回 False。

is_land(xpt, ypt)

xpt 和 ypt 是要进行计算的坐标点 1）坐标必须在地图坐标系中 2）在 Basemap 构造器中 resolution 不能是 None 3）如果点在 land 区域，将会使用 resolution 多边形进行计算，因此结果将依赖于 resolution
不能传入数组进行计算，只能一个点一个点的计算
使用 landpolygons 属性或 matplotlib PATH.contains_points方法进行多点计算 [注3]。上述方法不仅可以传递多个点，事实上，可以从 shp 文件中获取多个多边形进行判断。见 Basemap系列教程：使用shapefiles绘制地图 填充多边形 部分

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

map = Basemap(projection='aeqd', lon_0 = 10, lat_0 = 50, resolution='h')

x, y = map(0, 0)

print map.is_land(x, y)

#Idea taken from this post at StackOverflow: http://stackoverflow.com/questions/13796315/plot-only-on-continent-in-matplotlib/13811775#13811775
#I've re-written it to work with modern matplotlib versions

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from matplotlib.path import Path
import numpy as np

map = Basemap(projection='aeqd', lon_0 = 10, lat_0 = 50, resolution='h')

lons = [0., 0., 16., 76.]
lats = [0., 41., 19., 51.]

x, y = map(lons, lats)
locations = np.c_[x, y]

polygons = [Path(p.boundary) for p in map.landpolygons]

result = np.zeros(len(locations), dtype=bool) 

for polygon in polygons:
    result += np.array(polygon.contains_points(locations))

print result

location 是包含投影点的 numpy 数组
PATH 对象用于计算每一多边形
对于每一个PATH，使用 contains_points 计算每一个点。结果添加到 numpy 数组中，如果有一个多边形包含此点，结果将为 True

makegrid

makegrid 方法会创建随机格点，这些格点在地图坐标系中是等间距的点。用来获取使用地图投影构成的等间距格点的经度和纬度。

makegrid(nx, ny, returnxy=False)

nx 和 ny 定义的输出网格的大小。[注1]
如果 returnxy = True，地图投影中格点的位置将被返回
返回二维经纬度 numpy 数组，如果设置　returnxy = True，将返回地图坐标中的 x 和 y 数组

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

fig=plt.figure(figsize=(9, 3))

map = Basemap(width=12000000,height=8000000,
            resolution='l',projection='stere',
            lat_ts=50,lat_0=50,lon_0=-107.)

lons, lats, x, y = map.makegrid(30, 30, returnxy=True)

ax = fig.add_subplot(121)
ax.set_title('The regular grid')
map.scatter(x, y, marker='o')
map.drawcoastlines()


ax = fig.add_subplot(122)
ax.set_title('Projection changed')

map = Basemap(width=12000000,height=9000000,projection='aeqd',
            lat_0=50.,lon_0=-105.)
x, y = map(lons, lats)
map.scatter(x, y, marker='o')
map.drawcoastlines()

plt.show()

创建了两个地图。其中一个地图在规则网格中使用一种投影绘图，另一个地图使用另一种投影展示如何使用创建的经纬度矩阵
makegrid 函数中 returnxy 参数设置为 True 1）第一个地图中，直接使用 x 和 y 数组，因此具有相同的投影。如预期结果一致，创建了 30x30的规则网格 2）第二个地图使用的 lon 和 lat 矩阵，使用 basemap 实例重新投影数据到新的投影上。散点图显示了网格是不规则的

maskoceans

利用数据mask ocean 和 lake。

mpl_toolkits.basemap.maskoceans(lonsin, latsin, datain, inlands=True, resolution=’l’, grid=5)

此函数不是 basemap 实例方法，而是basemap模块中的独立函数
lonin 和 latin 是点的位置（二维数组）。此方法只能使用 latlon 投影。
datain 是包含要进行 mask 的数据数组
inland 控制是否也要 mask 内陆湖（lake）
resolution 设置要进行mask 的边界线精度，然后进行mask。默认为 'l'，与basemap默认值不同
grid 控制网格精度。可选值为10，5，2.5，1.25
下例展示了 resolution 和 grid 参数的差异

from mpl_toolkits.basemap import Basemap
from mpl_toolkits.basemap import maskoceans
from mpl_toolkits.basemap import interp
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., 
              resolution="h")

ds = gdal.Open("../sample_files/wrf.tiff")
lons = ds.GetRasterBand(4).ReadAsArray()
lats = ds.GetRasterBand(5).ReadAsArray()
data = ds.GetRasterBand(3).ReadAsArray()

x, y = map(lons, lats)

plt.figure(0)

mdata = maskoceans(lons, lats, data)

map.drawcoastlines(color = '0.15')
map.contourf(x, y, mdata)

plt.figure(1)

x2 = np.linspace(x[0][0],x[0][-1],x.shape[1]*5)
y2 = np.linspace(y[0][0],y[-1][0],y.shape[0]*5)

x2, y2 = np.meshgrid(x2, y2)

data2 = interp(data,  x[0], np.flipud(y[:, 0]), x2, np.flipud(y2),order=1)

lons2, lats2 = map(x2, y2, inverse=True)
mdata = maskoceans(lons2, lats2, data2, resolution = 'c', grid = 10, inlands=True)

map.drawcoastlines(color = '0.15')
map.contourf(x2, y2, mdata)

plt.figure(2)

mdata = maskoceans(lons2, lats2, data2, resolution = 'h', grid = 10, inlands=True)

map.drawcoastlines(color = '0.15')
map.contourf(x2, y2, mdata)

plt.figure(3)

mdata = maskoceans(lons2, lats2, data2, resolution = 'h', grid = 1.25, inlands=True)

map.drawcoastlines(color = '0.15')
map.contourf(x2, y2, mdata)

plt.show()

第一个例子（第22行）直接进行 mask，图中可以看出结果比较粗糙，这是由于输入数据的精度低，而不是 maskoceans 参数设置的原因。

第二个例子显示的图更精细一些（29-36行），从而避免了由于数据导致的大像素块。可查看 interp 部分获取更多细节。但 maskoceans 函数仍使用的 grid 和 resolution 参数的最低配置。

注意 Florida 的 lake 并没有被 mask，这是因为数据粗糙的原因导致。

在此例中，resolution 设置为 'h'，lakes已经被 mask了，但 Florida 海岸仍然显示一些像素块（就是有锯齿存在）。

最后，当 grid 参数也设置为最高精度后， Florida 海岸也能显示的很好（译注：其实放大后看，仍存在锯齿）。

nightshade

绘制指定日期地图中的夜晚地区。

nightshade(date, color=’k’, delta=0.25, alpha=0.5, ax=None, zorder=2)

date 是 datetime.datetime 对象
color 设置阴影区域颜色
dalta 是阴影区精度。默认为 0.25，当时用更小值时很容易出问题
alpha 是透明度设置
zorder 可以改变图层位置

下例显示了 van der Grinten 投影的阴影区 [注2]

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

map = Basemap(projection='vandg',lon_0=0,resolution='c')

map.drawmapboundary(fill_color="#7777ff")
map.fillcontinents(color="#ddaa66",lake_color="#7777ff")

map.drawcoastlines()

map.nightshade(datetime.now(), delta=0.2)

plt.show()

注1：原文为 nx and n define the size of the output grid

注2： http://matplotlib.org/basemap/users/vandg.html

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2017-07-26，如有侵权请联系 cloudcommunity@tencent.com 删除

编程算法

numpy

本文分享自气象杂货铺微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体分享计划，欢迎热爱写作的你一起参与！

编程算法

numpy

登录后参与评论

0 条评论

热度