OpenStreetMap

OpenStreetMap（简称OSM）是一个网上地图协作计划，目标是创造一个内容自由且能让所有人编辑的世界地图。由于大众参与的持续修正以及采用了更加可靠的数据源，OSM的海岸线精度和分辨率比GMT开源岸线GSHHG更高。

German FOSSGIS（https://www.fossgis.de/） 已经从OSM海岸线中制作了shapefile格式的文件，下载地址为https://osmdata.openstreetmap.de/data/land-polygons.html。 注意下载WGS-84投影的Large polygons not split版本，文件大小约600M。下载后可以利用开源GIS软件QGIS软件快速打开数据集，并导出感兴趣的区域为文本格式，从而可以使用别的绘图软件进行绘图（如GMT）。

使用QGIS提取区域shp

下面是基本的数据选取和导出步骤。

• 读入全球的shp文件
• 缩放到目标区域
• toggle edit
• Edit-Split feature,绘制一条线分割出目标区域。
• Select feature
• save selected feature as shp

如果是导出某一个闭合的多边形，如海岛，则不用分割，直接选择select feature，然后在layers中右键选择export-save selected feature as

如果想要导出一般海岛（或任意区域）的房屋、道路等信息，可以从https://www.openstreetmap.org 网站操作，缩放到感兴趣区域，然后点击左上角的Export，导出osm格式文件，使用QGIS打开后，可以转换为shp格式。

shp转文本格式或者kml

保证已经安装GMT，然后命令行：

1
2

ogr2ogr -f OGR_GMT sy.gmt sy.shp
gmt gmt2kml sy.gmt  -Wthick,white -Fl >sy.kml

使用GMT导出GSHHG海岸线

GMT可以绘制地图，也可以导出岸线数据，并基于岸线进行空间分析。下面是简单的岸线导出、绘图和转换代码。

1
2
3
4
5
6

ps=example_border.ps
gmt psbasemap  -JM4.5i -K -Bag `gmt gmtinfo sy.gmt -I0.000001`> $ps
gmt pscoast -M -W1 -Df -R >cn.txt # -Df means full resolution
gmt psxy -R -J cn.txt -W0.4p,red -O -K>> $ps
gmt psxy -R -J sy.gmt -W0.8p,black -O>> $ps
gmt gmt2kml cn.txt  -Wthick,white -Fl >cn.kml

效果

白线为GSHHG岸线，红线为OSM岸线。表明OSM更为准确，在一些对岸线精度较高的场景中可以使用OSM。

其他参考：
GMT提供了一个例子：https://docs.generic-mapping-tools.org/dev/gallery/ex51.html#example-51
QGIS：https://www.qgis.org/en/site/

PROJ 转ITRF地心坐标到经纬度

1
2

$ echo 4097216.6805  4429119.0287 -2065771.3676 0 | cct -I +proj=cart +ellps=GRS80
 47.2292116930  -19.0183060232     1552.9924        0.0000

PROJ 转经纬度到ITRF地心坐标

1
2

$ echo 47.2292116930  -19.0183060232     1552.9924 0 | cct +proj=cart +ellps=GRS80
4097216.6805   4429119.0287  -2065771.3676        0.0000

文件的批量转换

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

 $ cat points_xyz.dat | head
2919786.0        -5383745.0      1774604.0
4097216.6805     4429119.0287    -2065771.3676
6347491.3770     -22944.8191     622823.1488
4913652.948      3945922.493     995383.145
+2765120.9       -4449250.25     -3626405.6
-3530185.489     4118797.337     3344036.931
4696990.000      723994.000      4239678.000
-2341332.8840    -3539049.5090   4745791.3600
918129.40        -4346071.20     4561977.80
-4052051.767     4212836.215     -2545106.027

 $ cat points_xyz.dat | cct -I +proj=cart +ellps=GRS80 | head
-61.5275339092   16.2622989639      -25.6724           inf
 47.2292116930  -19.0183060232     1552.9924           inf
 -0.2071110447    5.6414771807       83.4568           inf
 38.7663016596    9.0351340011     2439.1542           inf
-58.1398674816  -34.8737126913       42.0708           inf
130.5995927060   31.8240613970      314.6402           inf
  8.7626079453   41.9274508625       98.7752           inf
-123.4874696602   48.3897817140       31.7442           inf
-78.0713675911   45.9558003990      200.8294           inf
133.8855132190  -23.6701238732      603.3581           inf

注意到第四列是时间信息，如果文件中缺少，则转换结果中提示inf。

卫星测高通用工具GAT

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

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

计划加入的功能：

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

这里列举一些常用的，又不容易记忆的linux工具命令，以及window快捷键

当前目录打开cmd

1
2
3
4
5

在地址栏中直接输入 cmd

    在文件资源管理器中导航到你希望打开 CMD 的目录。
    在窗口的地址栏中，输入 cmd 然后按 Enter。
    命令提示符将会在该目录中打开。

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调用的一个例子源码：