Шаг 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: [-53.0004 -53 -52.5017 -52.7963 -52.0434 -52.0838 ... ]
                 Lon: [-73.3617 -73.3626 -72.8909 -73.7034 -73.7392 ... ]
               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.shp');
shapewrite(L1, shapepath)

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

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

southAmerica = readgeotable(shapepath,'CoordinateSystemType','geographic')

southAmerica=79×13 table
       Shape         South      North      West      East        Area       Level    LevelString    NumPoints    FormatVersi    Source    CrossesGree    GSHHS_ID
    ____________    _______    _______    ______    ______    __________    _____    ___________    _________    ___________    ______    ___________    ________

    geopolyshape      -53.9     71.994    191.89    325.21    3.7652e+07      1        "land"          971            3         "WVS"          0             1   
    geopolyshape    -55.055    -52.469    287.98    294.89         46269      1        "land"           19            3         "WVS"          0            42   
    geopolyshape     -52.36    -51.235    300.26    302.29        4402.4      1        "land"           11            3         "WVS"          0            85   
    geopolyshape    -55.675    -54.933    289.99    292.02        2669.1      1        "land"           10            3         "WVS"          0            97   
    geopolyshape    -50.065     -48.67    284.53    285.62        4021.9      1        "land"            9            3         "WVS"          0           115   
    geopolyshape     -52.26    -51.336    298.92    300.81        4326.2      1        "land"            9            3         "WVS"          0           121   
    geopolyshape    -54.135     -53.38    286.33    287.85        3809.7      1        "land"            8            3         "WVS"          0           143   
    geopolyshape    -43.448    -41.768    285.55    286.65         10188      1        "land"            8            3         "WVS"          0           144   
    geopolyshape    -53.478     -52.73       287    288.46        3006.8      1        "land"            7            3         "WVS"          0           151   
    geopolyshape     10.041     10.841    298.07    299.09          5174      1        "land"            7            3         "WVS"          0           174   
    geopolyshape    -54.276    -53.566     289.1     289.8        1420.3      1        "land"            5            3         "WVS"          0           292   
    geopolyshape    -44.937    -44.364    286.54    287.31         955.5      1        "land"            5            3         "WVS"          0           329   
    geopolyshape    -54.401    -53.936    288.29    289.04         704.5      1        "land"            5            3         "WVS"          0           330   
    geopolyshape    -51.125     -50.71    285.04    285.64           462      1        "land"            5            3         "WVS"          0           339   
    geopolyshape    -50.458    -50.001    284.54    285.24         612.9      1        "land"            5            3         "WVS"          0           341   
    geopolyshape    -48.711    -48.137    284.98    285.53        1039.8      1        "land"            5            3         "WVS"          0           343   
      ⋮

Обратите внимание на то, что два самых долгих имен полей, '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

Ссылка

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

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

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

Кредиты

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

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

  >> type gshhs_c.txt