Документация

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

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

workingFolder = tempdir;

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

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

Используйте функцию gunzip MATLAB, чтобы распаковать 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 для Южной Америки

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

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

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

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 включает четыре уровня береговых линий:

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

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

ans = 1×3

     1     2     3

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

S(104).LevelString

ans = 
'lake'

показывает, что показывают 104, озеро (функция Уровня 2).

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

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

Этот пример управляет лучшими двумя уровнями иерархии GSHHS, вставляя каждое "озеро" в массу прилегающей земли.

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

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, если пара ограничительных рамок пересекается и ложь в противном случае). Входные параметры 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: сохраните результаты в файле форм

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

gshhs_c_SouthAmerica.shp
gshhs_c_SouthAmerica.shx
gshhs_c_SouthAmerica.dbf

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

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

Шаг 8: подтвердите файл форм

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

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

southAmerica = 79x1 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, формат) раньше хранил атрибуты в файлах форм. (В целом при записи файлов форм можно хотеть переопределить имена полей дольше, чем 11 символов во избежание или управлять эффектами автоматического усечения.)

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

delete([shapepath '.*'])

Отобразите геоstruct, импортированный из нового файла форм. Отметьте различные "дыры" в озерах указания полигона Южной Америки и береговых линиях других расширенных масс воды во внутренней части континента.

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

Ссылка

Вессел, P. и В. Х. Ф. Смит, 1996, глобальная последовательная, иерархическая, база данных береговой линии с высоким разрешением, Журнал Геофизического Исследования, Издания 101, стр 8741-8743.

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

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

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

Кредиты

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

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

  >> type gshhs_c.txt