gravity-toolkit安装和使用笔记整理一

安装

软件

安装anaconda和jupyter notebook，并在conda中创建虚拟环境，命名为grace，创建虚拟环境见[[Anaconda安装与环境配置]]

工具包

https://github.com/tsutterley/gravity-toolkit

https://github.com/tsutterley/geoid-toolkit/

https://github.com/tsutterley/read-GRACE-geocenter

环境

ubuntu16.04

依赖

conda激活grace环境，安装依赖。

激活：conda activate grace

安装：conda install -c conda-forge package_name

根目录(gravity-toolkit)

参照github作者中的gravity-toolkit/README.rst#dependencies安装依赖，包括
参照github作者中的gravity-toolkit/requirements.txt安装第三方库，包括
- boto3
- future
- lxml
- matplotlib
- numpy
- python-dateutil
- pyyaml
- scipy

.binder文件夹

github中作者中的gravity-toolkit/.binder/environment.yml中的依赖，详细信息如下。

name: gravity_toolkit
channels:

conda-forge

dependencies:

boto3
cartopy
future
h5py
ipywidgets>=7.6,<8.0(=7.7.5)
ipympl
jupyter
jupyterlab=3
jupyterlab_widgets
lxml
markupsafe==2.0.1
matplotlib
netCDF4
notebook
numpy
python>=3.6
python-dateutil
pyyaml
scipy
tk
pip
pip:
- git+https://github.com/tsutterley/geoid-toolkit.git

.github/workflows

github中作者中的gravity-toolkit/.github/workflows/python-request.yml中的安装依赖，详细信息如下。

此工作流将安装 Python 依赖项、运行测试并使用各种 Python 版本进行 lint。

name: Python on pull request

on:
pull_request:

jobs:
build:

runs-on: ${{ matrix.os }}
strategy:
  matrix:
    os: [ubuntu-latest, macos-latest, windows-latest]
    python-version: [3.8]
env:
  OS: ${{ matrix.os }}
  PYTHON: ${{ matrix.python-version }}
defaults:
  run:
    shell: bash -l {0}

steps:
- uses: actions/checkout@v2
- name: Set up conda ${{ matrix.python-version }}
  uses: conda-incubator/setup-miniconda@v2
  with:
    auto-update-conda: true
    python-version: ${{ matrix.python-version }}
    activate-environment: gravity_toolkit
    environment-file: .binder/environment.yml
- name: Create conda Test Environment
  run: |
    conda install flake8 pytest pytest-cov cython
    pip install git+https://github.com/tsutterley/read-GRACE-geocenter.git
- name: Lint with flake8
  run: |
    # stop the build if there are Python syntax errors or undefined names
    flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics
    # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide
    flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
- name: Test with pytest
  run: |
    pip install --no-deps .
    pytest --verbose --capture=no --cov=./ --cov-report=xml \
        --username=${{ secrets.EARTHDATA_USERNAME }} \
        --password=${{ secrets.EARTHDATA_PASSWORD }} \
        --webdav=${{ secrets.PODAAC_PASSWORD }}
- name: Archive code coverage results
  uses: actions/upload-artifact@v2
  with:
    name: code-coverage-report
    path: ./coverage.xml

.doc文件夹

github中作者中的gravity-toolkit/doc/environment.yml安装doc 环境，包括下列库。安装命令：conda install -c conda-forge package_name

name: gravity-docs
channels:
- conda-forge
  dependencies:
- docutils<0.18
- fontconfig
- fonttools
- freetype
- future
- graphviz
- lxml
- markupsafe==2.0.1
- numpy
- numpydoc
- pip
- python>=3.6
- python-dateutil
- pyyaml
- scipy
- sphinx
- sphinx_rtd_theme
- pip:
  - sphinx-argparse>=0.4
  - ..
github中作者中的gravity-toolkit/doc/source/getting_started/Citations.rst中的依赖，包括下列库。

安装

pip 安装(或选择下面的setup安装)

python3 -m pip install –user git+https://github.com/tsutterley/read-GRACE-geocenter.git

python3 -m pip install –user git+https://github.com/tsutterley/geoid-toolkit.git

python3 -m pip install –user git+https://github.com/tsutterley/gravity-toolkit

setup安装

进入gravity-toolkit文件夹（setup.py）同级目录，运行python setup.py install

同理，将read-GRACE-geocenter和geoid-toolkit拷贝入gravity-toolkit文件夹，运行python setup.py install

拷贝

将github上下载的文件夹中的script文件夹拷贝到~/anaconda3/envs/grace/lib/python3.8/site-packages/gravity_toolkit-1.2.0-py3.8.egg/gravity_toolkit 文件夹下便于直接访问script下的函数。

环境变量设置

将gravity_toolkit、gravity_toolkit/scripts和geoid_toolkit添加到环境变量(~~/.bashrc)~~中。

# <<< GRACE conda initialize <<<
alias grace="source activate grace"
export PATH="~/anaconda3/envs/grace/lib/python3.8/site-packages/gravity_toolkit-1.2.0-py3.8.egg/gravity_toolkit:$PATH"
export PATH="~/anaconda3/envs/grace/lib/python3.8/site-packages/gravity_toolkit-1.2.0-py3.8.egg/gravity_toolkit/scripts:$PATH"
export PATH="~/anaconda3/envs/grace/lib/python3.8/site-packages/geoid_toolkit-1.1.1-py3.8.egg/geoid_toolkit:$PATH"

添加完成后，终端执行source ~/.bashrc

测试

终端执行cnes_grace_sync.py --help和podaac_cumulus.py --help

1	podaac_cumulus.py --directory="/home/manasakura/draft" --release RL06

备注

setup安装后，在conda虚拟环境中就已存在該工具包，路径为~/anaconda3/envs/grace/lib/python3.8/site-packages，在调用工具包中的函数时实际上是调用site-packages中的函数，而不是从github上下载的文件夹中的函数，所以如果按照github文件夹的路径进行函数访问将不会访问到。

学习笔记

文档学习

doc文件夹

doc学习顺序可按照该文档进行

函数用法可查看readdocs.io和api_reference

gravity_toolkit文件夹

SLR文件夹

将C20.py, CS2.py, C30.py, C40.py, C50.py计算出的球协系数值（去除AOD）与matlab中计算出的进行对比，二者答案一致（只校对了GRACE_months 1～4月，对应2002年1月～2002年4月），matlab的计算值保存为~/coding/gravity_workspace/myworkspace/cs2m_c30_c40_c50_gracemonth_1_4_check.mat

C20.py

C20_RL06.txt(PROC=CSR)

SLR.C20中的data变量不是直接读取的C20_RL06.txt中的Column 2中的C20，而是将第二列的C20数据减去Column 5中的AOD产品，详细信息查看原文件中的comment如下。

The RL-06 background gravity model includes solid earth and ocean tides, solid earth and ocean pole tides, and the Atmosphere-Ocean De-aliasing (AOD1B) product. However, the monthly mean of the AOD model has been restored; subtract Column 5 to remove it.

源代码如下：

dinput['time'] = file_input['time']
dinput['month'] = gravity_toolkit.time.calendar_to_grace(dinput['time'])
# monthly spherical harmonic replacement solutions
dinput['data'] = file_input['C20'].copy()
# monthly spherical harmonic formal standard deviations
dinput['error'] = file_input['sigma']*1e-10
# Background gravity model includes solid earth and ocean tides, solid
# earth and ocean pole tides, and the Atmosphere-Ocean De-aliasing
# product. The monthly mean of the AOD model has been restored.
if AOD:
    # Removing AOD product that was restored in the solution
    dinput['data'] -= file_input['AOD']*1e-10

CS2.py

C21(22)_S21(22)_RL06.txt(PROC=CSR)`

读取的CS21/CS22是文件中SLR CS21/CS22（Column 2,3）减去AOD(大气-海洋去混叠,Column 6,7）之后的CS21/CS22,AOD简述可查看Atmosphere-Ocean De-aliasing For GRACE-FO - NASA/ADS 和 Evaluation of high-frequency oceanographic signal in GRACE data: Implications for de-aliasing

大气-海洋去锯齿 (AOD) 产品是一组 Stokes 系数，表示由于大气-海洋质量分布变化引起的非潮汐引力的时间效应。它被用作处理数据流的力模型，以防止高频大气和海洋运动的影响混叠到月解中。

源代码如下：

dinput['time'] = content['time'].copy()
dinput['month'] = gravity_toolkit.time.calendar_to_grace(dinput['time'])
# remove the monthly mean of the AOD model
dinput['C2m'] = content['C2'] - content['C2aod']*10**-10
dinput['S2m'] = content['S2'] - content['S2aod']*10**-10
# scale SLR solution sigmas
dinput['eC2m'] = content['eC2']*10**-10
dinput['eS2m'] = content['eS2']*10**-10

C30.py

CSR_Monthly_5x5_Gravity_Harmonics.txt(PROC=CSR)

利用read_SLR_harmonics.py计算C30完整球协系数，并返回C30完整球协系数，GRACE_month为从2002.0起算的月份（见time.py中的calendar_to_grace函数）。从GGMO5C平均重力场模型（头文件中）和SLR观测值（正文，系数以相对于平均重力场的变化形式给出，已经去除AOD）相加求C30完整球协系数，并将求得的C30完整球协系数作为返回值。

核心代码如下：

#### CSR_Monthly_5x5_Gravity_Harmonics.txt
# GGM05C重力场表头（头文件）
    L    M         C                  S               sigma C      sigma S
# SLR观测值（正文内容，球协系数异常值）表头及单位
n, m, Cnm, Snm, Cnm+AOD, Snm+AOD, C-sigma, S-sigma, MJD_mid_point, Year_mid_point
Units 1.0e-10
#### read_GRACE_harmonics.py
# fill mean field Ylms（头文件中GGMO5C平均重力场模型）
mean_Ylms['clm'][l1,m1] = np.float64(line[2].replace('D','E'))
mean_Ylms['slm'][l1,m1] = np.float64(line[3].replace('D','E'))
mean_Ylm_error['clm'][l1,m1] = np.float64(line[4].replace('D','E'))
mean_Ylm_error['slm'][l1,m1] = np.float64(line[5].replace('D','E'))

# fill anomaly field Ylms and rescale to output（SLR观测的cs_lm异常值和cs_lm_error）
Ylm_anomalies['clm'][l1,m1,d] = np.float64(line[2])*SCALE
Ylm_anomalies['slm'][l1,m1,d] = np.float64(line[3])*SCALE
Ylm_anomaly_error['clm'][l1,m1,d] = np.float64(line[6])*SCALE
Ylm_anomaly_error['slm'][l1,m1,d] = np.float64(line[7])*SCALE

# calculate full coefficients and full errors（C30完整球协系数）
Ylms['clm'][:,:,d] = Ylm_anomalies['clm'][:,:,d] + mean_Ylms['clm'][:,:]
Ylms['slm'][:,:,d] = Ylm_anomalies['slm'][:,:,d] + mean_Ylms['slm'][:,:]
Ylms['error']['clm'][:,:,d]=np.sqrt(Ylm_anomaly_error['clm'][:,:,d]**2 +
    mean_Ylm_error['clm'][:,:]**2)
Ylms['error']['slm'][:,:,d]=np.sqrt(Ylm_anomaly_error['slm'][:,:,d]**2 +
    mean_Ylm_error['slm'][:,:]**2)

#### C30.py
# CSR 5x5 + 6,1 file from CSR and extract C3,0 coefficients
Ylms = gravity_toolkit.read_SLR_harmonics(SLR_file, HEADER=True)
# extract dates, C30 harmonics and errors
dinput['time'] = Ylms['time'].copy()
dinput['data'] = Ylms['clm'][3,0,:].copy()
dinput['error'] = Ylms['error']['clm'][3,0,:].copy()

C40.py

CSR_Monthly_5x5_Gravity_Harmonics.txt(PROC=CSR)

同C30.py的方法计算C40完整球协系数。

#### CSR_Monthly_5x5_Gravity_Harmonics.txt 同上C30.py

#### read_GRACE_harmonics.py 同上C30.py

#### C30.py
# read 5x5 + 6,1 file from CSR and extract C4,0 coefficients
Ylms = gravity_toolkit.read_SLR_harmonics(SLR_file, HEADER=True)
# extract dates, C40 harmonics and errors
dinput['time'] = Ylms['time'].copy()
dinput['data'] = Ylms['clm'][4,0,:].copy()
dinput['error'] = Ylms['error']['clm'][4,0,:].copy()

C50.py

CSR_Monthly_5x5_Gravity_Harmonics.txt(PROC=CSR)

同C30.py的方法计算C50完整球协系数。

#### CSR_Monthly_5x5_Gravity_Harmonics.txt 同上C30.py

#### read_GRACE_harmonics.py 同上C30.py#### C50.py

# read 5x5 + 6,1 file from GSFC and extract coefficients
Ylms = gravity_toolkit.read_SLR_harmonics(SLR_file, HEADER=True)
# calculate 28-day moving-average solution from 7-day arcs
dinput.update(gravity_toolkit.convert_weekly(Ylms['time'],Ylms['clm'][5,0,:], DATE=DATE, NEIGHBORS=28))
# no estimated spherical harmonic errors
dinput['error'] = np.zeros_like(DATE,dtype='f8')

gen_averaging_kernel.py

生成平均核系数，从而最大限度地减少总误差。

其中输入参数为

INPUTS:
gclm: cosine spherical harmonics of exact averaging kernel

gslm: sine spherical harmonics of exact averaging kernel

eclm: measurement error in the cosine harmonics

eslm: measurement error in the sine harmonics

sigma: variance of the surface mass signal

hw: Gaussian radius of the kernel in kilometers

需要用到的计算公式如下：

区域核函数（exact averaging kernel）

$$
\vartheta(\theta, \phi)=0 \quad outside \quad the \quad basin\
\vartheta(\theta, \phi)=1 \quad inside \quad the \quad basin
$$

求gclm/gslm

$$
\begin{Bmatrix}
\Delta{\vartheta}{lm}^c \ \Delta\hat{\vartheta}{lm}^s
\end{Bmatrix}=
\int \vartheta(\theta,\phi) \widetilde{P}_{lm}(cos\theta)
\begin{Bmatrix}
cosm\phi \ sinm\phi
\end{Bmatrix}d\Omega
$$

离散形式为

$$
\begin{Bmatrix}
\Delta{\vartheta}{lm}^c \ \Delta\hat{\vartheta}{lm}^s
\end{Bmatrix}=
\sum_{i=1}^{180} \sum_{j=1}^{360} \vartheta(\theta_i,\lambda_j)\widetilde P_{lm}(cos\theta_i)sin(m \lambda_j)

\begin{Bmatrix}
cosm\phi \ sinm\phi
\end{Bmatrix}
sin \theta_i \Delta \theta_i \Delta \lambda_j
$$

附上$\vartheta(\theta, \phi)$的球谐展开形式

$$
\vartheta(\theta,\phi) = \frac{1}{4\pi}

\sum_{l=0}^{\infty}\sum_{m=0}^l \widetilde{P}{lm}(cos\theta)(\vartheta{lm}^c \Delta C+\vartheta_{lm}^s \Delta S)
$$

计算Wlms：

gen_disc_load.py

计算均匀 disc load（负荷，与mass含义相近）的重力（gravitational）球谐系数。

gen_harmonics.py

gen_harmonics.py函数进行球谐系数展开主要是利用下列公式：

$$
\Delta C_{lm}=\frac{1}{4\pi} \sum_{i=1}^{180} \sum_{j=1}^{360}\Delta \sigma(\theta_i,\lambda_j)P_{lm}(cos\theta_i)cos(m \lambda_j)sin \theta_i \Delta \theta_i \Delta \lambda_j
$$

$$
\Delta S_{lm}=\frac{1}{4\pi} \sum_{i=1}^{180} \sum_{j=1}^{360}\Delta \sigma(\theta_i,\lambda_j)P_{lm}(cos\theta_i)sin(m \lambda_j)sin \theta_i \Delta \theta_i \Delta \lambda_j
$$

其中，

$$
\Delta \theta = \Delta \lambda = \pi/180
$$

上面公式与Chen（2015）中公式有所不同，Chen（2015）公式如下：

$$
\Delta C_{lm}=\frac{3}{4\pi \rho_E a} \frac{k_l + 1}{2l+1} \sum_{i=1}^{180} \sum_{j=1}^{360}\Delta \sigma(\theta_i,\lambda_j)P_{lm}(cos\theta_i)cos(m \lambda_j)sin \theta_i \Delta \theta_i \Delta \lambda_j
$$

$$
\Delta S_{lm}=\frac{3}{4\pi \rho_E a} \frac{k_l + 1}{2l+1} \sum_{i=1}^{180} \sum_{j=1}^{360}\Delta \sigma(\theta_i,\lambda_j)P_{lm}(cos\theta_i)sin(m \lambda_j)sin \theta_i \Delta \theta_i \Delta \lambda_j
$$

将空间域的数据转化为球协系数.

def demo2():
    grid_dir = '~/coding/gravity_workspace/gwrks0628/CSR_RL06_GSM_ICE6G-D_High_Res_wCSR_C20_wCSR_21_wCSR_22_wCSR_C30_wCSR_C40_wCSR_C50_wTellus_DEG1_cmwe_L60M60_r300km_FL_004-248.nc'
    grid = spatial().from_netCDF4(grid_dir)
    grid_data = grid.data[:, :, 0]
    lon = grid.lon
    lat = grid.lat
    Ylms = gen_harmonics(grid_data, lon, lat, METHOD='integration')
    print("Finished!")

在上述代码中，grid_dir路径指向的文件是由[GRACE-Spatial-Maps.ipynb]文件生成的grid对象，包含计算出的等效水高格网。gen_stokes函数可将該格网（二维）展开为球谐系数，其计算结果与matlab包中的grid2cs.m计算结果相同。matlab代码如下：

grid_dir = '~/coding/gravity_workspace/gwrks0628/CSR_RL06_GSM_ICE6G-D_High_Res_wCSR_C20_wCSR_21_wCSR_22_wCSR_C30_wCSR_C40_wCSR_C50_wTellus_DEG1_cmwe_L60M60_r300km_FL_004-248.nc';;
info = ncinfo(grid_dir);
variables = info.Variables;
grid_data = ncread(grid_dir, 'z');    % n_lon * n_lat * n_month
lon = ncread(grid_dir, 'lon');        % 西经（-180-0）：东经（0-180）
lat = ncread(grid_dir, 'lat');        % 北纬（90-0）：南纬（0：-90）
time = ncread(grid_dir, 'time');
% matlab包要求grid_data格网数据呈东经：西经（0：180：-180/0）、北纬：南纬（90：-90）分布。
clear grid_data_temp;
grid_data_temp(:,:) = grid_data(1,:,:);
grid_data_temp2 = grid_data_temp';
grid_data_temp3 = [grid_data_temp2(:,721:1440) grid_data_temp2(:,1:720)];        % n_lat * n_lon, 呈东经：西经（0：180：-180/0）、北纬：南纬（90：-90）分布
cs2 = gmt_grid2cs(grid_data_temp3/100, 60);

需要注意的是，在利用matlab包进行计算时，grid单位换算成了m（grid_data_temp3/100），上面两个计算结果相同，结果如下。

备注：该函数不涉及units类的引入和factor的计算。

gen_stokes.py

将空间域的数据转化为球协系数，可在ocean_stokes.py中查看其如何使用。

公式引用自Chen（2015），计算公式如下，

因为

$$
M_e = \rho_e V = \frac{4 \pi a^3 \rho_e}{3}
$$

上面公式等价于下式，

$$
\Delta C_{lm}=\frac{a^2}{M_e} \frac{k_l + 1}{2l+1} \sum_{i=1}^{180} \sum_{j=1}^{360}\Delta \sigma(\theta_i,\lambda_j)P_{lm}(cos\theta_i)cos(m \lambda_j)sin \theta_i \Delta \theta_i \Delta \lambda_j
$$

$$
\Delta S_{lm}=\frac{a^2}{M_e} \frac{k_l + 1}{2l+1} \sum_{i=1}^{180} \sum_{j=1}^{360}\Delta \sigma(\theta_i,\lambda_j)P_{lm}(cos\theta_i)sin(m \lambda_j)sin \theta_i \Delta \theta_i \Delta \lambda_j
$$

核心实现代码如下

    #### gen_stokes.py

    # longitude degree spacing in radians
    dphi = dlon*np.pi/180.0
    # colatitude degree spacing in radians
    dth = dlat*np.pi/180.0
    factors = gravity_toolkit.units(lmax=LMAX)

    # reforming data to lonXlat if input latXlon
    sz = np.shape(data)
    data = data.T if (sz[0] == nlat) else np.copy(data)

    elif (UNITS == 1):
    # Default Parameter: Input in cm w.e. (g/cm^2)
    dfactor = factors.spatial(*LOVE).cmwe
    int_fact[:] = np.sin(th)*dphi*dth

    # Calculating cos/sin of phi arrays
# output [m,phi]
m = np.arange(MMAX+1)
ccos = np.cos(np.dot(m[:,np.newaxis],phi))
ssin = np.sin(np.dot(m[:,np.newaxis],phi))

# Calculating fully-normalized Legendre Polynomials
# Output is plm[l,m,th]
plm = np.zeros((LMAX+1, MMAX+1, nlat))
# added option to precompute plms to improve computational speed
if PLM is None:
    # if plms are not pre-computed: calculate Legendre polynomials
    PLM, dPLM = plm_holmes(LMAX, np.cos(th))

# Multiplying by integration factors [sin(theta)*dtheta*dphi]
# truncate legendre polynomials to spherical harmonic order MMAX
for j in range(0,nlat):
    plm[:,m,j] = PLM[:,m,j]*int_fact[j]

# Initializing preliminary spherical harmonic matrices
yclm = np.zeros((LMAX+1, MMAX+1))
yslm = np.zeros((LMAX+1, MMAX+1))
# Initializing output spherical harmonic matrices
Ylms = gravity_toolkit.harmonics(lmax=LMAX, mmax=MMAX)
Ylms.clm = np.zeros((LMAX+1, MMAX+1))
Ylms.slm = np.zeros((LMAX+1, MMAX+1))
# Multiplying gridded data with sin/cos of m#phis
# This will sum through all phis in the dot product
# output [m,theta]
dcos = np.dot(ccos,data)
dsin = np.dot(ssin,data)
for l in range(LMIN,LMAX+1):# equivalent to LMIN:LMAX
    mm = np.min([MMAX,l])# truncate to MMAX if specified (if l > MMAX)
    m = np.arange(0,mm+1)# mm+1 elements between 0 and mm
    # Summing product of plms and data over all latitudes
    # axis=1 signifies the direction of the summation
    yclm[l,m] = np.sum(plm[l,m,:]*dcos[m,:], axis=1)
    yslm[l,m] = np.sum(plm[l,m,:]*dsin[m,:], axis=1)
    # Multiplying by factors to convert to fully normalized coefficients
    Ylms.clm[l,m] = dfactor[l]*yclm[l,m]
    Ylms.slm[l,m] = dfactor[l]*yslm[l,m]

# return the output spherical harmonics object
return Ylms

#### units.py
        # set default keyword arguments
        kwargs.setdefault('include_elastic', True)
        fraction = np.ones((self.lmax+1))
        # compensate for elastic deformation within the solid earth
        if kwargs['include_elastic']:
            fraction += kl[self.l]
        # degree dependent coefficients
        # norm, fully normalized spherical harmonics
        self.norm = np.ones((self.lmax+1))
        # cmwe, centimeters water equivalent [g/cm^2]
        self.cmwe = self.rho_e*self.rad_e*(2.0*self.l+1.0)/fraction/3.0
        # mmwe, millimeters water equivalent [kg/m^2]
        self.mmwe = 10.0*self.rho_e*self.rad_e*(2.0*self.l+1.0)/fraction/3.0

使用 gen_harmonics.py中使用过的同种grid格网数据，进行球谐系数展开，代码如下。

def load_lovenumbers():
    """
    read love numbers
    Returns
    -------
    hl, kl, ll: tuple
    """
    love_numbers_file = utilities.get_data_path(['data', 'love_numbers'])
    header = 2
    columns = ['l', 'hl', 'kl', 'll']
    hl, kl, ll = read_love_numbers(love_numbers_file, LMAX=60,
                                   HEADER=header, COLUMNS=columns, REFERENCE='CF', FORMAT='tuple')
    return hl, kl, ll

def demo():
    grid_dir = '~/coding/gravity_workspace/gwrks0628/CSR_RL06_GSM_ICE6G-D_High_Res_wCSR_C20_wCSR_21_wCSR_22_wCSR_C30_wCSR_C40_wCSR_C50_wTellus_DEG1_cmwe_L60M60_r300km_FL_004-248.nc'
    love_number = load_lovenumbers()
    grid = spatial().from_netCDF4(grid_dir)

    grid_data = grid.data[:,:,0]
    lon = grid.lon
    lat = grid.lat
    Ylms = gen_stokes(grid_data, lon, lat, LOVE=love_number)
    print("Finished!")

运行结果如下

grace_date.py

相关介绍查看test_gravity_notebook.ipynb。

grace_find_months.py

相关介绍查看test_gravity_notebook.ipynb。

grace_months_index.py

相关介绍查看test_gravity_notebook.ipynb。

grace_input_months.py

读取GRACE/GRACE- fo文件中指定的球谐l（阶数）和m（次数）以及指定的日期范围。

改正包括：

一阶项替换(如果指定)。

用单反值替换C20(如果指定)。

用179+月份的单反值替换C21/S21/C22/S22/C30/C50(如果指定)。

更正使用GAE, GAF和GAG文件的ECMWF大气“跳跃”。

Wahr et al.(2015)对极潮漂移的修正。

疑问：該函数中CS21/CS22/C30/C40/C50项为什么要设置当grace_month大于179（源码中为 grace_month > 176）时才能用SLR观测值替换GRACE观测值，请问是否在哪篇文献中有所提及，是否可以解答或者提供一下文献链接呢？代码如下（此处仅展示C21项的替换代码，S21/CS22/C30/C40/C50代码基本相同）。

for i,grace_month in enumerate(months):
    count = np.count_nonzero(C21_input['month'] == grace_month)
    if (count != 0) and (grace_month > 176):
        k, = np.nonzero(C21_input['month'] == grace_month)
        grace_Ylms['clm'][2,1,i] = np.copy(C21_input['C2m'][k])
        grace_Ylms['slm'][2,1,i] = np.copy(C21_input['S2m'][k])
        grace_Ylms['eclm'][2,1,i] = np.copy(C21_input['eC2m'][k])
        grace_Ylms['eslm'][2,1,i] = np.copy(C21_input['eS2m'][k])

geocenter.py

from_tellus()

读取一阶项文件TN-13_GEOC_CSR_RL06.txt，不做处理，详细用法见test_gravity_notebook.ipynb.

harmonics.py

drift(self, t, epoch=2003.3)

将球协系数速度场按照时间积分计算漂移。

核心代码

# # gravity_test_substep.py
GRACE_Ylms = harmonics().from_dict(Ylms)
GIA_Ylms_rate = gia(lmax=LMAX).from_GIA(filename=GIA_file,
    GIA=GIA, mmax=MMAX)
GIA_Ylms = GIA_Ylms_rate.drift(GRACE_Ylms.time, epoch=2003.3)

# # drift(self, t, epoch=2003.3)
# calculate drift
for i,ti in enumerate(t):
    temp.clm[:,:,i] = self.clm*(ti - epoch)
    temp.slm[:,:,i] = self.slm*(ti - epoch)
# assign degree and order fields
temp.update_dimensions()
return temp

from_dict(self, d, **kwargs)

将字典转化为harmonics类型。

用法

Ylms = grace_input_months(base_dir, PROC=PROC, DREL=DREL, DSET=DREL, LMAX=60, start_mon=4, end_mon=24,
                                 missing=[6, 7, 18], SLR_C20=SLR_C20, DEG1=DEG1, SLR_21=SLR_21, SLR_22=SLR_22,
                                 SLR_C30=SLR_C30, SLR_C40=SLR_C40, SLR_C50=SLR_C50)
GRACE_Ylms = harmonics().from_dict(Ylms)

mean(self, apply=False, indices=Ellipsis)

执行后就会更新GRACE_Ylms内的变量，返回值为球协系数均值。

用法

1	GRACE_Ylms.mean(apply=True)

ocean_stokes.py

读取land-sea mask文件并转化为球协系数。

ocean_stokes(LANDMASK, LMAX, MMAX=None, LOVE=None, VARNAME=’LSMASK’, SIMPLIFY=False)

读取landsea.nc文件

下载源：https://www.ncl.ucar.edu/Applications/Data/cdf/landsea.nc 。

1x1 degree mask distributed from UCAR as part of NCL。
landsea.nc文件中的变量名包括lon,lat和LSMASK三类，可使用ncview landsea.nc进行查看。
landsea.nc中各值的含义：
0： ocean
1：land
2：lake
3：small island
4：ice shelf

ocean_stokes代码使用如下：

from gravity_toolkit import ocean_stokes
landsea_mask_file = "/home/manasakura/coding/gravity_workspace/Data/LANDSEA_MASK/landsea.nc"
from gravity_toolkit.utilities import get_data_path
from gravity_toolkit.read_love_numbers import read_love_numbers
love_numbers_file = get_data_path(['data','love_numbers'])
hl,kl,ll = read_love_numbers(love_numbers_file, HEADER=2, FORMAT='tuple', REFERENCE='CF')
ocean_Ylms = ocean_stokes(LANDMASK=landsea_mask_file, LMAX=60, MMAX=60, LOVE=(hl,kl,ll), VARNAME='LSMASK',
    SIMPLIFY=False)
print(ocean_Ylms)

ocean_stokes函数会调用spatial.py中的spatial.from_netCDF4()函数，主要用于读取landsea.nc文件，代码如下：

1
2
3

# # ocean_stokes.py
landsea = spatial().from_netCDF4(LANDMASK,
        date=False, varname=VARNAME)

ocean_stokes函数还会调用gen_stokes.py中的gen_stokes()函数，用于将空间域的数据转化为球协系数，代码如下：

# # ocean_stokes.py
# create land function
nth,nphi = landsea.shape
land_function = np.zeros((nth,nphi), dtype=np.float64)
# combine land and island levels for land function
indx,indy = np.nonzero((landsea.data >= 1) & (landsea.data <= 3))
land_function[indx,indy] = 1.0
# remove isolated points if specified
if SIMPLIFY:
    land_function -= find_isolated_points(land_function)
# ocean function reciprocal of land function
ocean_function = 1.0 - land_function
# convert to spherical harmonics (1 cm w.e.)
Ylms = gen_stokes(ocean_function.T, landsea.lon, landsea.lat,
    UNITS=1, LMIN=0, LMAX=LMAX, MMAX=MMAX, LOVE=LOVE)

上一段代码可简单表达为制作海洋掩膜（海洋部分的值为1），并转化为球协系数。

read_GRACE_harmonics.py

读取GRACE/GRACE-FO Level-2 数据。用法见test_gravity_notebook.ipynb。

read_love_numbers.py

load love_numbers

Load Love Numbers match expected for reference frame。用法见test_gravity_notebook.ipynb。

load Load_Love2_CE.dat

Gegout (2005) Load Love Numbers。用法见test_gravity_notebook.ipynb。

load PREM-LLNs-truncated.dat

Wang et al. (2012) Load Love Numbers。用法见test_gravity_notebook.ipynb。

read_GIA_model.py

读取GIA文件（球协系数）。

read_GIA_model(input_file, GIA=None, MMAX=None, DATAFORM=None, **kwargs)

从文件中（ICE-6G_High_Res_Stokes_trend.txt）直接读取GIA球协系数值（GRACE Approximation）。用法见test_gravity_notebook.ipynb。

gia.from_GIA(self, filename, **kwargs)

从文件中（ICE-6G_High_Res_Stokes_trend.txt）直接读取GIA球协系数值（GRACE Approximation）。用法见test_gravity_notebook.ipynb。

spatial.py

对空间数据进行读、写等操作。

spatial.from_netCDF4(self, filename, **kwargs)

读取nc文件。

landsea.nc

下载：https://www.ncl.ucar.edu/Applications/Data/cdf/landsea.nc

landsea.nc文件中包含lon, lat, LSMASK三个变量，该文件可利用ncview可视化查看，查看命令ncview landsea.nc。

核心代码：

# # ocean_stokes.py
# VARNAME = "LSMASK"
# LANDMASK = "/home/manasakura/coding/gravity_workspace/Data/LANDSEA_MASK/landsea.nc"
landsea = spatial().from_netCDF4(LANDMASK,
        date=False, varname=VARNAME)

from_netCDF4()函数返回值为landsea（spatial类型），包含attributes, data, extent, filename, lon, lat, spacing等成员变量。data变量中包含的即是landsea.nc文件中的”LSMASK”变量值，下面需要对data变量多加以说明。

data：180×360的ndarray类型的矩阵，上半矩阵为南半球，下半矩阵为北半球，左半矩阵为东半球，右半矩阵为西半球。即是以0度经线和南纬90度纬线的交叉点为坐标原点建立坐标系，横坐标为i，纵坐标为j，i越往东值越大，j越往北值越大，用numpy读取該矩阵也是同理，按照此坐标系可确定矩阵每个坐标（i，j）所对应的地理坐标。

scripts文件夹

time.py

calendar_to_grace(year,month=1,around=np.floor)

**comments:**convert calendar dates to GRACE/GRACE-FO months

1	grace_month = around(12.0*(year - 2002.0)) + month

#GRACE #toolkit

gravity-toolkit安装和使用笔记整理一

https://singyutang.github.io/2023/06/03/gravity-toolkit安装和使用笔记整理一/

作者

SingyuTang

发布于

2023年6月3日

许可协议

测试文章上一篇

利用Clash科学上网下一篇