Skip to content

Commit

Permalink
Add example of buffer (#1080)
Browse files Browse the repository at this point in the history
  • Loading branch information
ZMAlt authored Aug 4, 2023
1 parent 7f33e3f commit 876baad
Show file tree
Hide file tree
Showing 3 changed files with 104 additions and 0 deletions.
Binary file added source/examples/ex008/coast_buffer.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
43 changes: 43 additions & 0 deletions source/examples/ex008/ex008.sh
Original file line number Diff line number Diff line change
@@ -0,0 +1,43 @@
#!/usr/bin/env bash

cat > testline1.txt << EOF
>
7.8 50.4
7.8 50.6
8.1 50.6
8.1 50.4
7.8 50.4
EOF

cat > testline2.txt << EOF
>
8 50
9.5 50.5
EOF

cat > testline3.txt << EOF
>
9 49.6
9 50
9.8 50
9.8 49.8
9.4 49.8
9.4 49.6
EOF

gmt spatial testline1.txt -Sb0.05 > testline1_buffer.txt
gmt spatial testline2.txt -Sb0.05 > testline2_buffer.txt
gmt spatial testline3.txt -Sb0.05 > testline3_buffer.txt

gmt begin buffer png
gmt basemap -R7.5/10/49.5/51 -JX15c -B

gmt plot testline1.txt -W1p,green
gmt plot testline1_buffer.txt -W1p,red

gmt plot testline2.txt -W1p,green
gmt plot testline2_buffer.txt -W1p,red

gmt plot testline3.txt -W1p,green
gmt plot testline3_buffer.txt -W1p,red
gmt end show
61 changes: 61 additions & 0 deletions source/examples/ex008/index.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,61 @@
:orphan:

绘制缓冲区
===========

:示例贡献者: |周茂|

----

GMT 的 :doc:`/module/gmtspatial` 通过调用 `Geos 库 <https://libgeos.org/>`__ 可实现部分 GIS
分析功能,包括生成缓冲区。下面对笛卡尔坐标多边形/线数据生成缓冲区。

.. gmtplot:: ex008.sh
:width: 50%
:show-code: true

从上图看出,在某些生成的缓冲区中,存在一些异常点。经过分析,推测该 bug 生成原因为:GMT 调用 GEOS
生成缓冲区多边形后,又调用取最大外部环函数将其从多边形转换为线数据导致部分点丢失,丢失的部分点即
图中出现的异常点。目前该 bug 还未精确定位和修复。在实际使用中,如果生成缓冲区不是明显的凹多边形,且该
部分丢失的区域不影响显示或后续分析,用户可手动删除位于缓冲区文件尾部的异常点。更加稳健的做法为使用
GEOS 的 julia/python 库实现上述功能,上述封装无需多边形到线数据的转换,因此不存在该 bug。此外,还
可使用 GMT 的 julia 封装 `GMT.jl <https://github.com/GenericMappingTools/GMT.jl>`__ 中的 `buffer`
函数,该函数通过 GDAL 调用 GEOS,同样不存在上述 bug。

目前,缓冲区生成算法是基于笛卡尔坐标实现的,因此会将输入的地理坐标假定为笛卡尔
坐标,这在使用地理坐标时,通常是不合适的。正常的做法是,先将多边形坐标投影,然后计算
缓冲区,最后再将缓冲区做投影逆变换得到地理坐标,其中的投影变换或逆变换可通过
:doc:`/module/mapproject` 实现,选择的投影最好为等面积投影。下面分别展示不使用投影
和使用投影的地理坐标缓冲区计算。

.. gmtplot::
:caption: 海岸线缓冲区
:width: 60%

gmt begin coast_buffer

# 提取边界
gmt coast -ETW -M -Dc > tw.geo
# 适当删除面积太小的岛屿等特征,以防生成缓冲区计算太慢
gmt spatial -Qc1000+h tw.geo -fg > tw_temp.geo
# 生成缓冲区,宽度为 0.5 度
gmt spatial tw_temp.geo -Sb0.5 > tw_buffer.geo

# 坐标投影
gmt mapproject -R116/124/20/26 -Jh1:1 -C -F tw_temp.geo > tw_temp.car
# 生成缓冲区, 0.5 度约为 56 km
gmt spatial tw_temp.car -Sb56000 > tw_buffer.car
# 逆投影
gmt mapproject -R116/124/20/26 -Jh1:1 -C -F -I tw_buffer.car > tw_buffer_proj.geo

gmt plot -R116/124/20/26 -JQ5i tw.geo -W0.5p,black -B
gmt plot tw_buffer.geo -W0.5p,red
gmt plot tw_buffer_proj.geo -W0.5p,blue

gmt end show

通过手动编辑缓冲区文件,可得到下图。

.. figure:: coast_buffer.png
:width: 60%
:align: center

0 comments on commit 876baad

Please sign in to comment.