Шаг 1: Определение рабочей папки

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

workingFolder = tempdir;

Шаг 2: GNU ® разархивировать и индексировать уровень GSHHG с грубым разрешением

GSHHG доступен в широком диапазоне пространственных разрешений. В этом примере используются данные с наименьшим разрешением из двоичного файла gshhs_c.b. ZIP-копия этого файла содержится в папке 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» панели инструментов картографирования:

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, вставляя каждое «озеро» в окружающий земельный массив.

Экстракт GSHHS уровня 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 в массив 3-D:

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

Проверьте каждую пару элементов Level 1 - Level 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 '.*'])

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

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, Глобальная самосогласованная, иерархическая база данных береговой линии высокого разрешения, Journal of Geophysical Research, Vol. 101, pp. 8741-8743.

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

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

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

Кредиты

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

Для получения дополнительных сведений выполните следующее:

  >> type gshhs_c.txt