Шаг 1: Задайте рабочую папку

Этот пример создает несколько временных файлов и использует переменную workingFolder для обозначения их местоположения. Используемое здесь значение определяется выходом tempdir команда, но вы можете легко настроить это.

workingFolder = tempdir;

Шаг 2: GNU ® Unzip и индекс слоя GSHHG грубого разрешения

GSHHG доступен в широкой области значений пространственных разрешений. Этот пример использует данные с самым низким разрешением из двоичного файла gshhs_c.b. zipped копия GNU этого файла включена в папку Mapping Toolbox™ data в пути MATLAB.

Используйте MATLAB gunzip функция для декомпрессии gshhs_c.b.gz и создайте файл gshhs_c.b в местоположении, обозначенном workingFolder. Затем создайте файл индекса, gshhs_c.i, в той же папке. В целом, наличие файла индекса помогает ускорить последующие вызовы gshhs функция. Обратите внимание, что при использовании 'createindex' опция, gshhs не извлекает данные.

files = gunzip('gshhs_c.b.gz', workingFolder);
filename = files{1};
indexfile = gshhs(filename, 'createindex');

Шаг 3: Импорт данных GSHHG для Южной Америки

Выберите данные для определенного четырехугольника широты-долготы и импортируйте его как массив «geostruct» Mapping Toolbox:

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

Если вы закончили извлечение данных, можно удалить декомпрессированный файл GSHS и файл индекса.

delete(filename)
delete(indexfile)

Шаг 4: Исследуйте набор данных

Исследуйте первый элемент массива геомассивов структур S. В дополнение к Lat и Lon массивы координат, обратите внимание на различные поля атрибутов, которые присутствуют.

S(1)

ans = struct with fields:
            Geometry: 'Polygon'
         BoundingBox: [2x2 double]
                 Lat: [1x972 double]
                 Lon: [1x972 double]
               South: -53.9004
               North: 71.9942
                West: 191.8947
                East: 325.2054
                Area: 3.7652e+07
               Level: 1
         LevelString: 'land'
           NumPoints: 971
       FormatVersion: 3
              Source: 'WVS'
    CrossesGreenwich: 0
            GSHHS_ID: 1

GSHHS включает четыре уровня береговой линии:

Проверьте, какие уровни включают импортированные данные. The Level поле содержит числовые номера уровней.

levels = [S.Level];
unique(levels)

ans = 1×3

     1     2     3

The LevelString поле обеспечивает их интерпретацию. Для примера,

S(104).LevelString

ans = 
'lake'

Показы этого признака 104 являются озером (функция уровня 2).

В этом примере из-за низкого разрешения или пространственного подмножества нет функций уровня 4.

Шаг 5: Извлечение двух верхних уровней в отдельные массивы геомассивов структур

Этот пример манипулирует двумя верхними уровнями иерархии GSHS, вставляя каждое «озеро» в окружающую наземную массу.

Извлечение GSHHS Level 1 (внешние береговые линии континентов и океанических островов):

L1 = S(levels == 1);

Извлечение уровня 2 (береговые линии озер и морей на многоугольниках уровня 1):

L2 = S(levels == 2);

Чтобы увидеть их пространственные зависимости, можно отобразить ребра уровня 1 как синие линии и ребра уровня 2 как красные линии:

figure
axesm('mercator', 'MapLatLimit', latlim, 'MapLonLimit', lonlim)
gridm; mlabel; plabel
geoshow([L1.Lat], [L1.Lon], 'Color', 'blue')
geoshow([L2.Lat], [L2.Lon], 'Color', 'red')
tightmap

Шаг 6: Объедините многоугольники уровня 2 в уровень 1

Задайте анонимную функцию предиката, чтобы обнаружить пересечения ограничивающих прямоугольников (возвращая true, если пара ограничивающих прямоугольников пересекается и false в противном случае). Входные параметры A и B являются ограничивающими прямоугольными матрицами вида 2 на 2

  [min(lon)  min(lat)
   max(lon)  max(lat)].

boxesIntersect = ...
    @(A,B) (~(any(A(2,:) < B(1,:)) || any(B(2,:) < A(1,:))));

Для удобства закольцовывания по ним скопируйте ограничивающие рамки уровня 1 в трехмерный массив:

L1boxes = reshape([L1.BoundingBox],[2 2 numel(L1)]);

Проверьте каждую пару функций уровня 1 - уровня 2 на предмет возможного пересечения. Посмотрите, polybool возвращает любой выход или нет, но избегает вызова polybool если сначала не обнаружено пересечение ограничивающего прямоугольника:

for k = 1:numel(L2)
    for j = 1:numel(L1)
        % See if bounding boxes intersect
        if boxesIntersect(L2(k).BoundingBox, L1boxes(:,:,j))
            % See if actual features intersect
            if ~isempty(polybool('intersection', ...
                L2(k).Lon, L2(k).Lat, L1(j).Lon, L1(j).Lat))
                % Reverse level 2 vertex order before merge to
                % correctly orient inner rings
                L1(j).Lon = [L1(j).Lon fliplr(L2(k).Lon) NaN];
                L1(j).Lat = [L1(j).Lat fliplr(L2(k).Lat) NaN];
            end
        end
    end
end

Шаг 7: Сохраните результаты в Shapefile

С одним вызовом на shapewrite, можно создать трио файлов,

gshhs_c_SouthAmerica.shp
gshhs_c_SouthAmerica.shx
gshhs_c_SouthAmerica.dbf

в рабочей папке.

shapepath = fullfile(workingFolder,'gshhs_c_SouthAmerica');
shapewrite(L1, shapepath)

Шаг 8: Проверьте Shapefile

Чтобы подтвердить результаты shapewrite, считайте новый файл shapefile в массив геомассивов структур southAmerica:

southAmerica = shaperead(shapepath, 'UseGeoCoords', true)

southAmerica=79×1 struct array with fields:
    Geometry
    BoundingBox
    Lon
    Lat
    South
    North
    West
    East
    Area
    Level
    LevelString
    NumPoints
    FormatVersi
    Source
    CrossesGree
    GSHHS_ID

Обратите внимание, что два самых длинных имени поля, 'FormatVersion' и 'CrossesGreenwich', были усечены до 11 символов. Это произошло во время вызова на shapewrite и неизбежно из-за жесткого 11-символьного предела в таблицах xBASE (формат .DBF), используемого для хранения атрибутов в shapefiles. (В целом при записи shapefiles можно переопределить имена полей длиннее 11 символов, чтобы избежать или контролировать эффекты автоматического усечения.)

Вы можете удалить новые файлы shapefiles из рабочей папки. (Этот пример должен быть очищен после себя; в реальном приложении вы, вероятно, захотите опустить этот шаг.)

delete([shapepath '.*'])

Отображение геоstruct, импортированных из нового файла shapefile. Обратите внимание на различные « отверстия» в южноамериканском многоугольнике, указывающие озера и береговые линии других протяженных водоемов в глубине континента.

figure
ax = axesm('mercator', 'MapLatLimit', latlim, 'MapLonLimit', lonlim);
ax.Color = 'cyan';
gridm; mlabel; plabel
geoshow(southAmerica, 'FaceColor', [0.15 0.8 0.15])
tightmap

Ссылка

Wessel, P. and W. H. F. Smith, 1996, A Global self-consistent, ierarchical, high-resolution shoreline database, Journal of Geophysical Research, Vol. 101, pp. 8741-8743.

Дополнительные данные

Полный набор данных GSHHG может быть загружен с веб-сайта Национального управления океанических и атмосферных исследований (NOAA) США. Переходите по ссылкам из

https://www.mathworks.com/help/map/finding-geospatial-data.html

Кредиты

Файл данных GSHHG представлен в Mapping Toolbox от доктора Пола Весселя из Гавайского университета и доктора Уолтера Х. Ф. Смита из NOAA.

Для получения дополнительной информации выполните команду:

  >> type gshhs_c.txt