Python绘制地球的重力地图

Boule

boule中定义了多种参考椭球，可用于表示地球、火星等星体的重力分布。可通过pip安装

pip install boule

boule中已经定义的椭球如下

椭球	GRS80	WGS84	MARS	MERCURY	MOON	VENUS	VESTA
星体	地球	地球	火星	水星	月球	金星	灶神星

这些椭球可直接调用。下面以WGS84椭球为例，获取其纬度为45°，高50米处的重力

import boule as blg = bl.WGS84.normal_gravity(latitude=45, height=50)
print(f"{g:.2f} mGal")  # 980604.35 mGal

Gal是重力加速度单位，读作伽利略，1Gal=1cm/s${}^2$。上述代码的输出结果，表示在纬度为45°，海拔50米的地方，其重力加速度为9.8060435m/s${}^2$。

【normal_gravity】即为获取重力的函数，其输入参数latitude, height即纬度和高度，其支持数组输入，由此可得到水平面处，从南极点到北极点的重力场变化，如图所示

代码如下

import matplotlib.pyplot as plt
import numpy as npL = np.linspace(-90, 90, 181)
g = bl.WGS84.normal_gravity(latitude=L, height=0)/1000
plt.plot(L, g)
plt.xlabel("latitude/°")
plt.ylabel("gravity/Gal")
plt.show()

结合cartopy的绘图功能，即可得到地球表面任意位置的重力场。

重力地图

接下来绘制地球重力场，首先创建经纬度网格，然后将地球的重力场映射到经纬度网格中，效果如下。

绘制代码为

import verde as vd
import cartopy.crs as ccrsL, B = vd.grid_coordinates(region=[0, 360, -90, 90], spacing=0.5,)
g = bl.WGS84.normal_gravity(latitude=B, height=10_000)/1e5
grid = vd.make_xarray_grid((L, B), g, data_names="gravity")ax = plt.subplot(111,projection = ccrs.Mollweide())
plt.contourf(grid.easting, grid.northing, grid.gravity, 60,transform=ccrs.PlateCarree(), levels=100)
plt.colorbar()
plt.show()

【verde】是一个空间数据处理和插值模块，其在上述代码中，主要起到生成经纬度网格的作用。

【grid_coordinates】通过输入空间区域以及间隔，生成相应的经纬度数组L, B。然后根据此生成的纬度网格，生成海拔为1000m处的重力场。

【make_xarray_grid】用于生成网格，其内部的【easting】，【northing】分别为东向和北向的两组坐标。其内部的网格名称，则通过data_names来指定。

最后，通过contourf进行绘图，在绘图过程中，使用了坐标转换，将其解析为PlateCarree投影。

从ensaio下载重力数据

Ensaio在葡萄牙语中是随笔的意思，是一个用于下载开源数据集的python库。其底层基于Pooch来下载和管理数据。这些数据的预处理方法可以在github上找到：fatiando-data。

下面绘制Ensaio中提供的重力场数据，如下图所示

代码为

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import ensaio
import xarray as xr# 下载并加载全球重力场数据
fname = ensaio.fetch_earth_gravity(version=1)
data = xr.load_dataset(fname)# 提取重力数据
gravity = data["gravity"]/1e5# 创建绘图
ax = plt.subplot(111, projection=ccrs.Robinson())
ax.set_global()  # 设置为全球范围# 绘制重力数据
img = ax.contourf(gravity.longitude, gravity.latitude,gravity, transform=ccrs.PlateCarree(), levels=100)# 添加海岸线
ax.coastlines(resolution="110m", color="black", linewidth=0.8)
plt.colorbar(img, pad=0.05, aspect=50)
plt.show()