卫星测高通用工具GAT

我开源了卫星雷达高度计（卫星测高）的一些基本程序，取名The Generic Altimetry Tools，简称GAT。名字的选择特别向GMT致敬，GMT帮助了全世界无数的地球科学研究科研人员。
GAP程序面向卫星测高的科研人员，学生等群体，主要目的是对卫星测高数据进行数据处理，因此采用数据处理能力较强的matlab语言和PYTHON语言编写，并调用GMT的基本程序，现在可以实现如下功能：

• Netcdf文件格式的读取
• SSH，SLA的计算
• 卫星测高交叉点分析
• 卫星测高绝对定标
• 卫星测高的湿延迟评估
• 卫星测高河湖的水位计算
• 卫星测高波形读取和绘图

计划加入的功能：

• 卫星测高交叉点平差
• 卫星测高共线平差
• 卫星测高重力反演
• 卫星测高水深反演
• 卫星测高中尺度涡探测
• 卫星测高的海平面变化分析
• 卫星测高波形重定
• 卫星测高极地海冰厚度计算
• 卫星测高海冰识别
• 卫星测高极地波形分类

这里列举一些常用的，又不容易记忆的linux工具命令

zotero 添加SCIhub

1

Edit-Preference-Config editor-extensions.zotero.findPDFs.resolvers

添加

1
2
3
4
5
6
7
8
9

{
    "name":"Sci-Hub",
    "method":"GET",
    "url":"https://sci-hub.ren/{doi}",
    "mode":"html",
    "selector":"#pdf",
    "attribute":"src",
    "automatic":true
}

然后，使用Zotero Connector或者Find Available PDF均可下载已有的PDF文件。

nvidia

1
2
3
4
5
6
7
8

nvidia-smi # 查看显卡信息
watch -n 0.1 nvidia-smi #每隔0.1s更新
lspci|grep NVIDIA # 查看显卡驱动
lspci | grep -i vga
sudo sh NVIDIA-Linux-x86_64-440.64.run # 安装Nvidia驱动
# 或者下面代码安装cuda
sudo chmod +x cuda_11.0.3_450.51.06_linux.run
sudo ./cuda_11.0.3_450.51.06_linux.run

添加路径

1
2
3
4
5

sudo gedit ~/.bashrc
export PATH="/usr/local/cuda-11.0/bin:$PATH"
export LD_LIBRARY_PATH="/usr/local/cuda-11.0/lib64:$LD_LIBRARY_PATH"
source ~/.bashrc
nvcc -V # 查看cuda安装

删除cuda。（有时候需要删除后，重新安装别的版本）

1
2
3

cd /usr/local/cuda-11.0/bin/
sudo ./cuda-uninstaller
sudo rm -rf /usr/local/cuda-11.0

安装cupti

1
2

sudo apt-get install libcupti-dev
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/cuda/extras/CUPTI/lib64

Sphinx

1
2
3
4

conda install Sphinx
sphinx-build --version
sphinx-quickstart docs
sphinx-build -b html docs docs/_build/html

Colab

Colaboratory is built on top of Jupyter Notebook.

GPU可以快速提高训练的速度。Colab提供K80显卡，价格一万多元。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

# 查看显存情况
!/opt/bin/nvidia-smi

# 查看GPU
import tensorflow as tf
tf.test.gpu_device_name()
# 在使用哪个GPU
from tensorflow.python.client import device_lib
device_lib.list_local_devices()

# 查看CPU
!cat  /proc/cpuinfo

# RAM大小
!cat /proc/meminfo

使用google drive数据的方法：https://colab.research.google.com/notebooks/io.ipynb

1
2
3
4

from google.colab import drive
drive.mount('/content/drive')
import os
os.chdir('/content/drive/My Drive/')

Conda

Conda 比较慢，可以使用mamba加速：

1

conda install -c conda-forge mamba

之后，使用mamba替代conda使用。

Conda是一种python包管理工具，可以对环境进行封闭管理，不同环境之间不相互干扰。

1
2
3
4
5
6
7
8
9
10
11
12
13
14

conda --version # 版本号
conda env list # 列出所有环境
conda list -n python_3.9 # 列出非当前环境的包
conda list # 列出当前环境的包
conda create -n env_name package_name # 创建环境，如conda create -n python2 python=2.7 
conda activate python2.7 # 激活某环境
conda deactivate # 退出当前环境
conda install package_name #在当前环境中安装包
conda remove package #删除当前环境中的包
conda install --name env_name package_name #在指定环境中安装包
conda remove -- name env_name package 
conda create -n tensorflow_env tensorflow #安装tensorflow的CPU版本
conda install tensorflow # 在当前环境下安装tensorflow
conda search tensorflow # 查找版本

删除环境之前需要先deactivate

1
2

conda deactivate tensorflow_env
conda remove -n tensorflow_env --all

添加conda-forge channel并设置为最高优先级

1

conda config --add channels conda-forge

查看现在的channel状态和优先级

1

conda config --get channels

激活某个环境，在环境中设置channel的严格优先。注，如果该环境中已经安装了若干库，则需要先更新所有的库以保证大多数库从conda-forge下载，以保持库的一致性。

1
2
3
4
5
6

# activate my environment
conda activate my_env
# update all packages if needed
conda update --all
# set channel priority as strict
conda config --set channel_priority strict

使用添加频道 conda config –add 频道 —（你的频道） 您添加的最后一个频道获得最高优先级……所以请保持顺序。您可以添加频道，即使您已经拥有它们，以便更改优先级顺序.

恢复原来的源,删除添加的源和设置。

1

conda config --remove-key channels

Docker

理解image、container。image可以共用，比如Ubuntu环境，而且非常精简。container是每次启动后的封闭容器，一个Ubuntu可以启动无数个container，彼此不相干。

1
2
3
4
5
6
7
8
9
10

docker images
docker info
docker run -it ubuntu bash
docker ps
docker ps -aq
docker start (id)
docker restart
docker kill
docker kill $(docker ps -aq)
docker start -d (id)

退出容器：

1

exit

退出，但是还在后台运行，Ctrl+P+Q.
退出后再次进入docker start -d (id)或者docker exec -it 3e17550b10c8 bash。此外还有attach可以进入正在运行的窗口。

查看日志

1

docker logs -f -t --details 3e17550b10c8

hexo

下面是简单的上传博客文章的hexo命令，其中hexo d -m "message" 可以产生commit。

1
2
3

hexo new post "linux tools"
hexo g
hexo d -m "linux tools"

插入公式：

1
2

npm uninstall hexo-math --save
npm install hexo-renderer-mathjax --save

配置参考：

https://www.jianshu.com/p/9b9c241146bc

行内的公式用$公式内容$包裹起来，在公式中右键即可查看公式$\TeX$代码。
同样右键可查看公式代码。

例如：
\begin{equation} \sum_{k=1}^{n} \frac{1}{1+8 \sin ^{2} \frac{k \pi}{n}}=\frac{n\left(2^{n}+1\right)}{3\left(2^{n}-1\right)} \end{equation}

文章抬头需要加入类似：

1
2
3
4
5
6

---
title: hexo中插入数学公式
date: 2019-05-09 12:40:16
mathjax: true
tags: 
---

git

常用git命令

1
2
3
4

git status
git add .
git commit -m "add dac"
git push origin master

弹出openSSH的问题，要求第二次输入username和密码，注意这里需要登录GitHub，新建一个Personal access tokens。username就是note的名字，密码就是token。

更新代码到本地。如果本地没有修改，可以直接使用：

1

git pull

撤销不小心add的文件

1

git reset data/era5/4d/south/era5.pl.*.nc

.gitignore修改无效的问题
.gitignore只能忽略那些原来没有被track的文件，如果某些文件已经被纳入了版本管理中，则修改.gitignore是无效的，所以我们需要使用rm命令清除一下相关的缓存内容.这样文件将以未追踪的形式出现，对某个文件取消跟踪。

1
2
3

git rm -r --cached .
git add .
git commit -m 'update .gitignore'

git撤回commit

写完代码后，我们一般这样
git add . //添加所有文件
git commit -m “本功能全部完成”
执行完commit后，想撤回commit，怎么办？
这样凉拌：
git reset –soft HEAD^

不再追踪某文件

1

git rm --cached test/ja3_check/mypoints.kml

在VPN，比如clash，环境下使用，遇到错误：

1

fatal: unable to access 'https://github.com/yangleir/gnsswave.git/': failed to connect to github.com port 443: timed out

解决方案：

1
2

git config --global https.proxy http://127.0.0.1:7890
git config --global https.proxy https://127.0.0.1:7890

或者在用户根目录下，直接打开.gitconfig，手动修改：

1
2
3
4
5

[http]
	postBuffer = 524288000
	proxy = 127.0.0.1:7890
[https]
	proxy = 127.0.0.1:7890

批量解压gz

我下载了很多个文件夹，每个文件夹下都是gz压缩文件，现在想批量解压指定文件夹的文件。

文件下载如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

$ ls | head
cycle_0066/
cycle_0067/
cycle_0068/
cycle_0069/
cycle_0070/
cycle_0071/
cycle_0072/
cycle_0073/
cycle_0074/
cycle_0075/

$ ls cycle_0099 | head
global_sla_l2p_ntc_h2_C0099_P0001_20150704T020146_20150704T024914_20180712T075832.nc.gz
global_sla_l2p_ntc_h2_C0099_P0002_20150704T025439_20150704T034613_20180712T075832.nc.gz
global_sla_l2p_ntc_h2_C0099_P0003_20150704T034614_20150704T043339_20180712T075832.nc.gz
global_sla_l2p_ntc_h2_C0099_P0004_20150704T043901_20150704T053040_20180712T075832.nc.gz
global_sla_l2p_ntc_h2_C0099_P0005_20150704T053041_20150704T062248_20180712T075832.nc.gz
global_sla_l2p_ntc_h2_C0099_P0006_20150704T062321_20150704T071508_20180712T075833.nc.gz
global_sla_l2p_ntc_h2_C0099_P0007_20150704T071509_20150704T080720_20180712T075833.nc.gz
global_sla_l2p_ntc_h2_C0099_P0008_20150704T080722_20150704T085935_20180712T075833.nc.gz
global_sla_l2p_ntc_h2_C0099_P0009_20150704T085936_20150704T095149_20180712T075833.nc.gz
global_sla_l2p_ntc_h2_C0099_P0010_20150704T095150_20150704T104402_20180712T075833.nc.gz

解压cycle_0068和cycle_0069两个文件夹下的全部gz文件，并删除gz原始文件。

1

gzip -d cycle_006[8-9]/*.gz

批量转换图片格式

把pdf转为png，并设置density 300，表示图像每英寸面积内的像素点数，数值越高图片质量越高

1

for file in *.pdf; do magick convert -density 300 $file ${file%%.*}.png; done

视频编辑

截取视频

1
2

$ ffmpeg  -i ./linchang.mp4 -vcodec copy -acodec copy -ss 00:00:00 -to 00:00:16 ./cutout1.mp4 -y
$ ffmpeg  -i ./linchang.mp4 -vcodec copy -acodec copy -ss 00:00:29 -to 00:00:48 ./cutout2.mp4 -y

合并视频(对于MP4格式，需要先转ts再合并)

1
2
3

$ ffmpeg -i cutout1.mp4 -c copy -bsf:v h264_mp4toannexb -f mpegts input1.ts   
$ ffmpeg -i cutout2.mp4 -c copy -bsf:v h264_mp4toannexb -f mpegts input2.ts   
$ ffmpeg -i "concat:input1.ts|input2.ts" -c copy -bsf:a aac_adtstoasc -movflags +faststart output.mp4

慢速，三倍慢速

1

$ ffmpeg -i 2021-02-19-22-38-32.mp4 -vf  "setpts=3*PTS" test2.mp4

Latex转word

1

pandoc template.tex -o output.docx -w docx  --pdf-engine pdflatex

Linux 解压乱码问题

使用下面的编码即可。

1

unzip -O CP936 test.zip

上一篇使用了MATLAB模拟二维波动网格，使用GMT做谱分析和绘图。本文进一步简化模拟代码，仅使用GMT一种语言可完成波数谱的模拟，并展示一维和二维波数谱的差异。

GMT有几种计算波数谱的程序

GMT有spectrum1D和grdfft可以计算波数谱。

• spectrum1D:一维谱分析。
• grdfft:二维谱分析，含x,y和径向三个方向的计算。

模拟代码

下面利用完全GMT工具进行波数谱的模拟，首先利用grdmath模拟出一个规范的正弦波二维矩阵，形成网格文件。然后使用grdfft计算三个方向的谱，使用spectrum1D计算X方向的一维波数谱。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

#!/usr/bin/env bash
ps=testgrd.ps

# gmt grdmath -R-50/50/-50/50 -I0.2 X Y ADD SIN = x.nc
gmt grdmath -R-50/50/-50/50 -I0.2 X SIN = x.nc

gmt makecpt -Crainbow -T-1/1/0.1 > t.cpt
gmt grdimage x.nc -R-30/30/-30/30 -JX3i -Ct.cpt -Bxa -Bya -BWSne -P -K -Y7i > $ps

gmt grd2xyz x.nc | awk '{if ($2 == 0) print $0}' > x.txt
gmt psxy -R -J -O -K x.txt -Sc0.03i -Gblack >> $ps
# Do grdfft
gmt grdfft x.nc -Er+w+n -Gfft.txt
gmt grdfft x.nc -Ex+w+n -Gfftx.txt
gmt grdfft x.nc -Ey+w+n -Gffty.txt
gmt math fftx.txt -C1 SQR ffty.txt SQR ADD SQRT = fftrr.txt

# Plot specrum result from grdfft
gmt psbasemap -R1/20/1e-14/1e-2 -JX-2.8il/3il -Bxa1pf2g3+l"Wavelength" -Bya1pg+l"power spectrum" -BwNsE -K -X3.2i -O >> $ps
gmt psxy -R fft.txt  -J -W0.5p,red, -O -K >> $ps
gmt psxy -R fftx.txt  -J -W0.5p,blue, -O -K >> $ps
gmt psxy -R ffty.txt  -J -W0.5p,green, -O -K >> $ps
gmt psxy -R fftrr.txt  -J -W0.5p,black, -O -K >> $ps

gmt psbasemap -R0.05/1e0/1e-14/1e-2 -JX2.8il/3il -BSE -Bxa1p+l"wavenumber/cycle" -O -K>> $ps

gmt pslegend -D+w1.2i+jBL+o0.1i/0.1i -R -J -O -K -F+p1p+glightgray --FONT_ANNOT_PRIMARY=10p,1,black << EOF >> $ps
S 0.2i - 0.5i - 1p,red, 0.5i r
S 0.2i - 0.5i - 1p,blue 0.5i x
S 0.2i - 0.5i - 1p,green, 0.5i y
EOF

# Do spectrum1D
gmt psxy -R-30/30/-1/1 -JX3i/3i x.txt -i0,2 -Bxaf -Byaf -BWSne -O -K -X-3.2i -Y-4.5i >> $ps
gmt psxy -R -J x.txt -i0,2 -Sc0.03i -Gred -O -K >> $ps

awk '{print $1,$3}' x.txt | gmt spectrum1d  -S256 -W --GMT_FFT=brenner -N -i1 -D0.2  > fft_spec.txt
gmt psbasemap -R1/20/1e-12/1e2 -JX-2.8il/3il -Bxa1pf2g3+l"Wavelength" -Bya1pg+l"power spectrum" -BwNsE -K -X3.2i -O >> $ps
gmt psxy -R fft_spec.txt -J -W0.5p,red, -O -K >> $ps
gmt psbasemap -R0.05/1e0/1e-12/1e2 -JX2.8il/3il -BSE -Bxa1p+l"wavenumber/cycle" -O >> $ps

awk '{print $3}' x.txt| gmt gmtmath STDIN -Sl STD  SQR =

gmt psconvert  $ps -P -Tg -A

结果

z=sin(x)

上为二维FFT在三个方向的谱分析结果，黑色曲线为sqrt(x^2+y^2),径向的能量谱不与一维的平方和（开根号）相等。下为一维（彩图中横线的采样）谱分析结果（形状与二维FFT的X方向一致，但能量大小不一样）。

z=sin(x+y)

上图把函数改变，数据x和y方向周期依然是2pi,但是r方向是2pi*cos(pi/2)，即4.44。这样可以帮助理解r方向是什么方向。运行时，把下面的代码取消注释。

1

# gmt grdmath -R-50/50/-50/50 -I0.2 X Y ADD SIN = x.nc

注意的问题

• spectrum1D得到的波数谱符合能量守恒定律，即输入信号的方差等于功率谱密度的积分，也就是谱和X轴围起来的面积。例如正弦波的方差为0.5，而功率密度的积分也恰好为0.5。
• grdfft得到的波数谱不是密度谱，所以不遵循能量守恒。
• grdfft的三个方向的谱的关系是什么？r^2=x^2+y^2？

参考

wavenumber的模拟计算

阿基米德螺旋线

心形线

参考

参考

问：GMT自带的movie程序怎么制作动画？是否简单易用？

答：movie是GMT的一个动画制作程序，设计的比较复杂，上手时间长，不易用。希望尝试的话，可以运行自带例子，你会发现它设计的使用方法和经典GMT是差别很大的，个人学习的使用经历比较糟糕。

问：那怎么利用GMT做GIF动画？

答：可以使用ImageMagick工具，只需要一行命令即可。

1

magick convert -delay 64 -loop 0 *.png  out.gif

首先使用GMT制作每帧png，然后用上面的命令把连续的png图像制作成GIF动画。（GMT movie的基本原理可能也是如此，但设计的复杂度太高）

例子

这里以GNSS浮标观测的海面波动为例，希望展示波浪的时间变化规律。动画设置每一帧显示一个小时的波浪功率谱图像，不同帧使用不同的颜色表达。

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

#!/usr/bin/env bash
gmt set FORMAT_DATE_MAP "o dd" FORMAT_CLOCK_MAP hh:mm 
gnss1=./gnss_b1_hlf_allppk.txt
R=-R32400/119700/-0.565631/0.828242
# ==========================================================================================================
for((i=1;i<=23;i++));  
do   
	echo $i
	((name = $i + 100))
	ps=a00$name.ps
	r=$(($RANDOM%255))
	g=$(($RANDOM%255))
	b=$(($RANDOM%255))
	
	((max =  $i * 3600))
	((min = ($i-1) * 3600))
	echo $min

	gmt psbasemap -R2019-03-31T09:00:00/2019-04-01T09:00:00/-0.565631/0.828242 -Bpxa4Hf1h -Bsxa1Df1D -Bpya0.4g0.2f0.2+l"WSSE(HPF)/m" -JX10c/4c -BWSne -K > $ps
	gmt psbasemap -R2/20/1e-4/1e0 -JX-10cl/10cl -Bxa2f3g3+l"period:second" -Bya-1pg+l"power spectrum:m@+2@+/cycle/s" -BWNse -K  -Y2.7i -O >> $ps
	
	awk 'NR<'"$max"' && NR>'"$min"'  {print $1,$2}' $gnss1 | gmt spectrum1d  -S256 -W --GMT_FFT=brenner -N -i1 -D1 > pow5.txt
	awk ' NR<'"$max"'   {print $1,$2}' $gnss1| gmt psxy  $R  -JX10c/4c -W0.5p,lightgray -P -K -O -Y-2.7i >> $ps
	awk 'NR<'"$max"' && NR>'"$min"'  {print $1,$2}' $gnss1| gmt psxy  $R  -JX10c/4c -W0.5p,$r/$g/$b -P -K -O  >> $ps
	gmt psxy -R2/20/1e-4/1e0 -JX-10cl/10cl -K pow5.txt -W2.25p,$r/$g/$b -O -Y2.7i -i0,1 >> $ps

	gmt psbasemap -R0.05/0.5/1e-4/0.5 -JX10cl/10cl -Bxa2f3g3+l"Frequency:cycle/second" -BS  -O >> $ps
gmt psconvert $ps -A -Tg
done
magick convert -delay 64 -loop 0 *.png  out.gif

结果

在一张图上

如果是做PPT展示用，动画较好。如果是插入论文，则需要把不同时间的曲线叠加到一张图上。同样使用for循环，代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69

#!/usr/bin/env bash
# Testing gmt spectrum1d power spectrum values
gmt set FORMAT_DATE_MAP "o dd" FORMAT_CLOCK_MAP hh:mm 
ps=test.ps

# ==========================================================================================================
# Use GMT spectrum to calculate power spetrum of SSH
gnss1=./gnss_b1_hlf_allppk.txt
gnss2=./gnss_b2_hlf_allppk.txt

gmt psbasemap -R2019-03-31T09:00:00/2019-04-01T09:00:00/-0.565631/0.828242 -Bpxa4Hf1h -Bsxa1Df1D -Bpya0.4g0.2f0.2+l"WSSE(HPF)/m" -JX10c/4c -BWSne -K > $ps
gmt psbasemap -R2/20/1e-4/1e0 -JX-10cl/10cl -Bxa2f3g3+l"period:second" -Bya-1pg+l"power spectrum:m@+2@+/cycle/s" -BWNse -K  -Y2.7i -O >> $ps

R=-R32400/119700/-0.565631/0.828242
echo period > period1_wind.txt
# echo period > period_l.txt
# GNSS buoy1 set i<=23;
# GNSS buoy2 set i<=24;
for((i=1;i<=23;i++));  
do   
	echo $i
	r=$(($RANDOM%255))
	g=$(($RANDOM%255))
	b=$(($RANDOM%255))
	# ((a=1+2**3-4%3))
	((max =  $i * 3600))
	((min = ($i-1) * 3600))
	echo $min
	# echo $(expr $i \* 3 + 1);  
	# awk 中使用变量的小技巧“‘val’”
	
	awk 'NR<'"$max"' && NR>'"$min"'  {print $1,$2}' $gnss1 | gmt spectrum1d  -S256 -W --GMT_FFT=brenner -N -i1 -D1 > pow5.txt
	awk 'NR<'"$max"' && NR>'"$min"'  {print $1,$2}' $gnss1| gmt psxy  $R  -JX10c/4c -W0.5p,$r/$g/$b -P -K -O -Y-2.7i >> $ps
	gmt psxy -R2/20/1e-4/1e0 -JX-10cl/10cl -K pow5.txt -W2.25p,$r/$g/$b -O -Y2.7i -i0,1 >> $ps

	gmt gmtselect -R2/7/1e-4/1e0 pow5.txt | sort -k 2 | tail -n 1| awk '{print '"$i+8"',$1}' >> period1_wind.txt
done
gmt psbasemap -R0.05/0.5/1e-4/1e0 -JX10cl/10cl -Bxa2f3g3+l"Frequency:cycle/second" -BS  -O -K>> $ps

###B2
gmt psbasemap -R2019-03-31T09:00:00/2019-04-01T09:00:00/-0.565631/0.828242 -Bpxa4Hf1h -Bsxa1Df1D -Bpya0.4g0.2f0.2+l"WSSE(HPF)/m" -JX10c/4c -BwSnE -K  -X11c -Y-2.7i -O >> $ps
gmt psbasemap -R2/20/1e-4/1e0 -JX-10cl/10cl -Bxa2f3g3+l"period:second" -Bya-1pg+l"power spectrum:m@+2@+/cycle/s" -BwNsE -K  -Y2.7i -O >> $ps

R=-R32400/119700/-0.565631/0.828242
echo period > period2_wind.txt
# GNSS buoy1 set i<=23;
# GNSS buoy2 set i<=24;
for((i=1;i<=24;i++));  
do   
	echo $i
	r=$(($RANDOM%255))
	g=$(($RANDOM%255))
	b=$(($RANDOM%255))
	# ((a=1+2**3-4%3))
	((max =  $i * 3600))
	((min = ($i-1) * 3600))
	echo $min
	# echo $(expr $i \* 3 + 1);  
	# awk 中使用变量的小技巧“‘val’”
	
	awk 'NR<'"$max"' && NR>'"$min"'  {print $1,$2}' $gnss2 | gmt spectrum1d  -S256 -W --GMT_FFT=brenner -N -i1 -D1 > pow5.txt
	awk 'NR<'"$max"' && NR>'"$min"'  {print $1,$2}'  $gnss2| gmt psxy  $R  -JX10c/4c -W0.5p,$r/$g/$b -P -K -O -Y-2.7i >> $ps
	gmt psxy -R2/20/1e-4/1e0 -JX-10cl/10cl -K pow5.txt -W2.25p,$r/$g/$b -O -Y2.7i -i0,1 >> $ps
	gmt gmtselect -R2/7/1e-4/1e0 pow5.txt | sort -k 2 | tail -n 1| awk '{print '"$i+8"',$1}' >> period2_wind.txt
done
gmt psbasemap -R0.05/0.5/1e-4/1e0 -JX10cl/10cl -Bxa2f3g3+l"Frequency:cycle/second" -BS  -O >> $ps

gmt psconvert $ps -A -Tf
rm gmt.* h_f.txt smooth_track* t.d

获取数据

GNSS data 1
GNSS data 2

注意点

• GraphicsMagick是从ImageMagick中分离出来的，推荐下载ImageMagick。

问：我想计算海洋上的任意点到海岸线的（最短）距离，并制作一个网格文件存储这些距离，请问如何？

答：使用grdmath可获得全球任意点到海岸线距离的网格。这里的距离就是垂直最短距离。

问：我想绘制一条线，表示到海岸线的距离。

答：通过设置grdcontour可以导出指定的等值线坐标到文本，比如我们要绘制距离海岸线50km的线，或者专属经济区的线等。

问：我想根据到海岸线的距离提取点位，可以吗？

答：可以。下面的例子中有这个功能的使用。如从我国GNSS陆态网站点中提取出距离海岸线50km的点，并绘制。

在物理海洋学中有一个常用的量叫做wavenumber，翻译为波数，其频谱叫做波数谱。一般来说，遥感和测绘学科中对波数谱较为陌生，大家使用最多的时间相关的频率谱，而波数谱是空间谱。海洋学中定点观测（如浮标）获得的数据一般是时间序列，可以计算频率谱，而卫星（如高度计）得到的是沿轨数据，属于空间域，卫星数据可以用于计算波数谱。本文简单对此介绍，并使用GMT和MATLAB做简单的模拟，以求两种维度的谱概念清晰易懂。

wavenumber是什么？和frequency是什么关系？

• wavenumber是空间域频率，相当于单位长度有多少个波，也就是波数，其物理单位可以是CPKM(cycle per km)。 所谓波数谱分析就是指空间域的FFT分析，可以是一维的along-track分析，也可以是二维的空间分析，波数谱分析的Y轴单位一般是cm^2/CPKM，表示在某一波数上的能量。
• frequency是时间域的频率，相当于单位时间内的波的个数，大家都懂。
• wavelength是波长，就是一个波的空间长度（空间域）或者时间长度（时间域），时空通用。

  In the physical sciences, the wavenumber (also wave number or repetency) is the spatial frequency of a wave, measured in cycles per unit distance or radians per unit distance. Whereas temporal frequency can be thought of as the number of waves per unit time, wavenumber is the number of waves per unit distance.

• 维基百科wavenumber
• 参考图（来自维基百科）：
• k表示波数，是物理海洋中的一个常见符号。例如，波长为0.5km，那么线性的波数为2，表示单位长度（1km）内有两个波。k也可以用弧度表示,等于2Pi/wavelenth，线性波数2再乘以2Pi为4Pi，是以[rad/m]为单位的波数表达形式。

wavenumber怎么计算，MATLAB、PYTHON，GMT？

几乎所有的语言都可以计算波数（傅里叶分析），下列是几套较为成熟的开源软件：

• PYTHON: pyspec，论文中开源的代码和数据，含有详细的结果分析。
• Matlab: jlab，一套丰富的海洋大数据分析工具。
• GMT: GMT内置了两种计算PSD的程序，一维的spectrum1d和二维的grdfft。

MATALB+GMT的二维wavenumber模拟计算

MATLAB模拟二维网格

参数设置：

• 网格维度为1600*801，间隔为0.3m。
• 波长为50m和10m的叠加。
• 加入了白噪声。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57

% test NetCDF write
% Author: Yang Lei,2020-04-19,Qingdao. Email:leiyang@fio.org.cn
clear all;
delete ('test_simulate.nc')
lati = linspace(0,1600*0.3,1600);
long = linspace(0,1600*0.3,1600);

% simulate wave 
% build 1D sinewave
SineMid = sin(2*pi/10*long);

% translate 1D wave into 2D sinewave
SinewaveM = repmat(SineMid,[1600,1]);
imagesc(SinewaveM)
Angle = 10;
Irot = imrotate(SinewaveM, Angle, 'loose', 'bilinear');
% imagesc(Irot)

% figure out the image size of new Sinewave.
Start = floor(size(Irot,1)/4);
Stop = 3*Start;
 
% image crop
Icrop = Irot(Start:Stop,Start:Stop);
L_x=length(Icrop);

temp=Icrop;

ssh_wave = temp+0.1*rand(L_x,L_x);
imagesc(ssh_wave)

% netcdf.create
nccreate('test_simulate.nc','lat','Dimensions',{'lat' L_x});

ncwriteatt('test_simulate.nc', 'lat', 'standard_name', 'lat');
ncwriteatt('test_simulate.nc', 'lat', 'long_name', 'lat');
ncwriteatt('test_simulate.nc', 'lat', 'units', 'meter');
ncwriteatt('test_simulate.nc', 'lat', '_CoordinateAxisType', 'Lat');
% ncwriteatt('test_files.nc','/','standard_name','latitude');
nccreate('test_simulate.nc','lon','Dimensions',{'lon' L_x});
ncwriteatt('test_simulate.nc', 'lon', 'standard_name', 'lon');
ncwriteatt('test_simulate.nc', 'lon', 'long_name', 'lon');
ncwriteatt('test_simulate.nc', 'lon', 'units', 'meter');
ncwriteatt('test_simulate.nc', 'lon', '_CoordinateAxisType', 'Lon');

nccreate('test_simulate.nc','Q','datatype','double','Dimensions',{'lat' L_x 'lon' L_x});
ncwriteatt('test_simulate.nc', 'Q', 'standard_name', 'ssh');
ncwriteatt('test_simulate.nc', 'Q', 'long_name', 'sea surface height');
ncwriteatt('test_simulate.nc', 'Q', 'units', 'm');

% ncdisp('test_files.nc');
ncwrite('test_simulate.nc','lat',lati(1:L_x));
ncwrite('test_simulate.nc','lon',long(1:L_x));

ncwrite('test_simulate.nc','Q',ssh_wave);

ncdisp('test_simulate.nc');

利用GMT的grdfft进行二维FFT分析，得到波数谱。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

#!/bin/bash
gmt gmtset MAP_GRID_PEN_PRIMARY	= 0.1p,0/0/0,2_1_0.25_1:0

ps=gf-fft.ps

data_jizai=test_simulate.nc

gmt grdfft $data_jizai -Ex+w+n -Gx.txt
gmt grdfft $data_jizai -Ey+w+n -Gy.txt
gmt grdfft $data_jizai -Er+w+n -Gr.txt

awk '{print $1,$2}' x.txt > jizai.x 
awk '{print $1,$2}' y.txt > jizai.y
awk '{print $1,$2}' r.txt > jizai.z

x=`awk '{print $1,$2}' x.txt | sort -k 2 | head -n 1 | awk '{print $1}'`
y=`awk '{print $1,$2}' y.txt | sort -k 2 | head -n 1 | awk '{print $1}'`
r=`awk '{print $1,$2}' r.txt | sort -k 2 | head -n 1 | awk '{print $1}'`

echo $x $y $r
gmt gmtmath -Q $x $y DIV ATAN 180 pi DIV MUL 90 SUB NEG =

gmt psbasemap -R0.1/1e3/1e-15/1e5 -JX-4il/4il -Bxa1pf3g3+l"Wavelength/m" -Bya1pf1+l"power spectrum" -BWNe -K > $ps

gmt psxy -R jizai.x  -J -W0.5p,black, -O -K >> $ps
gmt psxy -R jizai.y -J -W0.5p,blue, -O -K >> $ps
gmt psxy -R jizai.z -J -W0.5p,red, -O -K >> $ps

gmt psbasemap -R1e-3/1e1/1e-15/1e5 -JX4il/4il -BSE -Bxa1p+l"wavenumber/cycle/m" -O -K>> $ps

gmt pslegend -D+w0.8i+jBL+o0i/.2i -R -J -O -F+p1p+glightgray --FONT_ANNOT_PRIMARY=7p << EOF >> $ps
S 0.2i - 0.5i - 1p,red 0.5i x
S 0.2i - 0.5i - 1p,blue 0.5i y  
S 0.2i - 0.5i - 1p,tan 0.5i r
EOF

gmt psconvert $ps -A -P -Tg

rm gmt.* jizai.* *.cpt *.d *.nc

结果计算得到三个方向的波长和径向传播角度：

1
2

61.478424015 10.2464040025 10.2464040025
9.46232220803

PLOT grid

1
2
3
4
5
6
7
8
9
10
11
12

#!/usr/bin/env bash
	
data=test_files.nc

ps=mss.ps
J=JX4i

gmt makecpt -Crainbow -T-13/13/0.01 -Z >mss.cpt
gmt psbasemap -R$data -$J -BWNSE -Bxya100 -K > $ps
gmt grdimage $data -$J -Cmss.cpt -K -O >> $ps
gmt psscale -DjCB+w1.5i+o0/-1.4c+h -Cmss -Bxaf -By+lm -R -J -O -V >> $ps
gmt psconvert mss.ps -A -P -Tg

Attention

• 如果matlab中的XY属性给的随意，写成了longitude，latitude，而实际上是笛卡尔坐标，则GMT在识别的时候仍旧认为是经纬度，因此不能大于360或者90。所以在写nc文件的时候要注意属性名称。
• kx和ky的计算分别使用-Ex和-Ey.
• 默认的计算方向是the radial direction
• 波的传播方向可以通过kx和ky的波长计算，atan(ky/kx)。
• 波的径向传播方向只能通过径向的FFT计算，不可以简单的相加kx和ky。

前期咨询了GMT开发者一些windows中使用GMT C API的问题，可以看这，发现直接使用GMT exe安装包中的库总是出现问题(可能和依赖有关系？和平台有关系？)，最终选择了使用visual studio 2019进行源码编译lib，并调用API试验。

API C调用的一个例子源码：

问答

问：上一篇文章尝试了使用cygwin来编译GMT，出现了一些问题。那么有没有别的办法编译GMT的github开发版（还没有正式发布的版本）？
答：根据官网指南，使用Visual Studio是一个可行的办法。

问：windows下编译源码复杂吗？需要依赖诸多库怎么办？
答：比Linux复杂。具体使用vcpkg来安装库。关于怎么安装和使用vcpkg，详见官网指南。vcpkg的安装较简单，使用的时候需要在线下载依赖库的压缩文件，程序自动下载的速度和网络有关，境内下载github源码较慢，这时可手动下载文件并复制到C:\vcpkg-master\downloads文件夹下，手动下载的地址即为命令窗口最下面的地址，然后根据C:\vcpkg-master\downloads\temp下面的文件名称对其重命名，再分次执行vcpkg install netcdf-c gdal pcre fftw3[core,threads] clapack openblas --triplet x64-windows，程序会根据自动分析是否安装和是否下载，安装好的库将跳过。

问答

问：在windows下通常都是直接安装GMT的exe程序，这样导致不能使用github上最新的开发版，而开发版通常fix了很多的bug。如果想在win上安装编译开发版本，怎么办？
答：应该尝试cygwin+cmake。

问：cygwin是什么？怎么使用？
答：可以看这里. cygwin的各个库升级都是通过setup**.exe实现的，每一次升级或者安装新的模块，都要打开这个程序。用习惯了，不觉得麻烦。

问：cmake怎么在这下面安装？
答：可以使用setup**.exe，选择cmake进行安装。

安装完成后打开命令行界面，查看cygwin版本：