Шаг 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=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, формат) раньше хранил атрибуты в файлах форм. (В общем случае при записи файлов форм можно хотеть переопределить имена полей дольше, чем 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