将海岸线数据 (GSHHG) 转换为 Shapefile 格式

“全球自洽分层高分辨率地理数据集” (GSHHG) [1] 是一套地理数据集，其中包含了海岸线、湖泊和主要河流的边界。GSHHG 数据可用于创建地图底图，以及进行需要地理搜索和几何量计算等操作的分析。

以下示例演示如何：

从 GSHHG 数据集提取一个子集。
使用多边形来表示数据。
将代表陆地和湖泊的多边形合并。
将合并后的数据保存为 shapefile 文件。

索引：低分辨率 GSHHG 数据

从一个 GNU 压缩文件中提取包含粗略 GSHHG 数据的文件。将文件作为地理数据结构数组读入工作区。

files = gunzip("gshhs_c.b.gz");
filename = files{1};

创建一个索引文件，以加快今后对 gshhs 函数的调用。请注意，当您使用 "createindex" 选项创建索引文件时，gshhs 函数不会提取数据。

indexfile = gshhs(filename,"createindex");

读取 GSHHG 数据

请指定南美洲周边区域的纬度和经度边界。然后，将该区域的 GSHHG 数据作为地理数据结构数组读入工作区。

latlim = [-60  15];
lonlim = [-90 -30];
S = gshhs(filename,latlim,lonlim);

使用 struct2geotable 函数将结构数组转换为地理空间表。地理空间表是一种 table 或 timetable 对象，其中包含一个 Shape 变量和属性变量。

GT = struct2geotable(S);

分析 GSHHG 数据

请注意，Shape 变量包含一个以地理坐标表示的多边形形状。

GT.Shape

ans=155×1 geopolyshape array with properties:
              NumRegions: [155×1 double]
                NumHoles: [155×1 double]
                Geometry: "polygon"
    CoordinateSystemType: "geographic"
           GeographicCRS: []
      ⋮

GSHHS 数据包含多个海岸线标高：

1 级多边形代表陆地区域。
第 2 级多边形代表湖泊。
第 3 级多边形代表湖中的岛屿。
第 4 级多边形表示位于湖泊中的岛屿内的池塘。

通过查询地理空间表中的 Level 变量，查找南美洲 GSHHG 数据中的海岸线高程。结果表明，南美洲 GSHHS 数据包含 1 级、2 级和 3 级多边形。

levels = GT.Level;
unique(levels)

此外，您还可以通过查询地理空间表中的 LevelString 变量，查阅南美洲 GSHHG 数据中海岸线高程的解释说明。该结果还表明，南美洲的 GSHHS 数据包含了湖中岛屿（第 3 级）、湖泊（第 2 级）和陆地（第 1 级）。

levelstrings = GT.LevelString;
unique(levelstrings)

ans = 3×1 cell
    {'island_in_lake'}
    {'lake'          }
    {'land'          }

显示陆地面积和湖泊

从地理空间表中创建两个子表。根据表示土地面积的表行，创建第一个子表。根据表中代表湖泊的行，创建第二个子表。

idxL1 = (levels == 1);
L1 = GT(idxL1,:);

idxL2 = (levels == 2);
L2 = GT(idxL2,:);

在地图上标出陆地面积和湖泊。将陆地区域用蓝色轮廓的多边形表示，将湖泊用红色轮廓的多边形表示。请注意，湖泊多边形还包括其他延伸的水域。

figure
worldmap(latlim,lonlim)
mlabel off
geoshow(L1,EdgeColor="#0072BD",FaceColor="none")
geoshow(L2,EdgeColor="#A2142F",FaceColor="none")

Figure contains an axes object. The hidden axes object contains 158 objects of type patch, line, text.

合并陆地区域和湖泊

将陆地区域多边形与湖泊多边形合并，使得合并后的结果是一组带有代表湖泊的孔洞的陆地区域多边形。

对于土地面积和湖泊地理空间表中的每一行：

获取陆地区域和湖泊多边形的边界框。
使用 boxesIntersect 辅助函数判断边界框是否重叠。此步骤通过省略对不重叠多边形形状的计算，从而提高了代码的性能。
如果边界框发生重叠，则通过计算陆地区域多边形与湖泊多边形的差集，在陆地区域多边形中生成孔洞。
在土地面积地理空间表中，将原始的土地面积多边形替换为更新后的土地面积多边形。

for landIDX = 1:height(L1)
    % land area bounding box
    landBBOX = L1.BoundingBox{landIDX};
    for lakeIDX = 1:height(L2)
        % lake bounding box
        lakeBBOX = L2.BoundingBox{lakeIDX};  
        % determine if bounding boxes overlap
        if boxesIntersect(landBBOX,lakeBBOX)
            landShape = L1.Shape(landIDX);
            lakeShape = L2.Shape(lakeIDX);
            % calculate difference of land area and lake polygons
            mergedShape = subtract(landShape,lakeShape);    
            % replace original land area polygon
            L1.Shape(landIDX) = mergedShape;                
        end
    end
end

将合并后的数据写入 Shapefile

将更新后的地理空间表写入 shapefile 文件。请注意，该 shapefile 由多个文件组成：一个主文件 (.shp)、一个索引文件 (.shx) 和一个属性文件 (.dbf)。

filename = "gshhs_c_SouthAmerica.shp";
shapewrite(L1,filename)

验证 Shapefile 文件

通过读取并显示文件中存储的数据，对 shapefile 进行验证。

将 shapefile 文件读入一个新的地理空间表中。

GT2 = readgeotable("gshhs_c_SouthAmerica.shp",CoordinateSystemType="geographic");

显示新表的变量名。请注意，由于属性文件 (.dbf) 规范中的字符限制，原始属性名 FormatVersion 和 CrossesGreenwich 被截断为 11 个字符。

GT2.Properties.VariableNames'

ans = 13×1 cell
    {'Shape'      }
    {'South'      }
    {'North'      }
    {'West'       }
    {'East'       }
    {'Area'       }
    {'Level'      }
    {'LevelString'}
    {'NumPoints'  }
    {'FormatVersi'}
    {'Source'     }
    {'CrossesGree'}
    {'GSHHS_ID'   }

将合并后的数据显示在地图上。请注意，陆地面积多边形中包含一些代表湖泊和其他大面积水体的空洞。

figure
worldmap(latlim,lonlim)
ax = gca;
setm(ax,"FFaceColor","#4DBEEE")
mlabel off
geoshow(GT2,FaceColor="#77AC30")

Figure contains an axes object. The hidden axes object contains 86 objects of type patch, line, text.

参考资料

[1] “GSHHG——一个全球性的、自洽的、分层的高分辨率地理数据库。”访问日期：2024 年 10 月 22 日。https://www.soest.hawaii.edu/pwessel/gshhg/。

辅助函数

辅助函数 boxesIntersect 返回逻辑值 1，表示由 bbox1 和 bbox2 指定的边界框相交；否则返回 0。bbox1 和 bbox2 是形式为 [lonmin latmin; lonmax latmax] 的 2×2 矩阵。

function tf = boxesIntersect(bbox1,bbox2)
    tf = ~(any(bbox1(2,:) < bbox2(1,:)) || any(bbox2(2,:) < bbox1(1,:)));
end