✍️ 周茂 • 📅 2023-04-04

# 绘制缓冲区

GMT 的 spatial 通过调用 Geos 库 可实现部分 GIS 分析功能，包括生成缓冲区。下面对笛卡尔坐标多边形/线数据生成缓冲区。

    #!/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

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

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

    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

海岸线缓冲区

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