文章目录
- 参考椭球
- 标准重力
- 重力地图
参考椭球
boule中定义了多种参考椭球，可用于表示地球、火星等星体的重力分布。可通过pip安装
pip install boule
安装完成后可直接调用
import boule as bl
boule中已经定义的椭球如下
| 椭球 | 星体 |
|---|---|
| GRS80 | 地球 |
| WGS84 | 地球 |
| MARS | 火星 |
| MERCURY | 水星 |
| MOON | 月球 |
| VENUS | 金星 |
| VESTA | 灶神星 |
这些椭球可直接调用
print(bl.WGS84)
Ellipsoid(name='WGS84', semimajor_axis=6378137, flattening=0.0033528106647474805, geocentric_grav_const=398600441800000.0, angular_velocity=7.292115e-05, long_name='World Geodetic System 1984', reference='Hofmann-Wellenhof, B., & Moritz, H. (2006). Physical Geodesy (2nd, corr. ed. 2006 edition ed.). Wien\u202f; New York: Springer.')
标准重力
下面以WGS84椭球为例，绘制地球表面的重力场分布。在Ellipsoid对象中，通过normal_gravity方法，可以获取指定纬度和高度处的重力，示例如下
import boule as bl
gamma = bl.WGS84.normal_gravity(latitude=45, height=50)
print(f"{gamma:.2f} mGal")
# 980604.35 mGal
Gal是重力加速度单位，读作伽利略，1Gal=1cm/s 2 ^2 2。上述代码的输出结果，表示在纬度为45°，海拔50米的地方，其重力加速度为9.8060435m/s 2 ^2 2。
其两个参数latitude和height可以数组，由此可得到不同纬度的重力加速度变化
import matplotlib.pyplot as plt
L = np.linspace(-90, 90, 181)
gamma = bl.WGS84.normal_gravity(latitude=L, height=0)
plt.plot(L, gamma)
plt.show()
重力地图
接下来绘制地球重力场，首先创建经纬度网格，然后将地球的重力场映射到经纬度网格中。
import numpy as np
import verde as vd
L, B = vd.grid_coordinates(region=[0, 360, -90, 90], spacing=0.5,)
gamma = bl.WGS84.normal_gravity(latitude=B, height=10_000)
grid = vd.make_xarray_grid((L, B), gamma, data_names="gravity")
然后调用pygmt中的绘图函数
import pygmt
fig = pygmt.Figure()
fig.grdimage(grid.gravity, projection="W20c", cmap="viridis")
fig.basemap(frame=["af", "WEsn"])
fig.colorbar(position="JCB+w10c", frame=["af", 'y+l"mGal"', 'x+l"WGS84"'])
fig.show()
