В этом примере показано, как к:
Извлеките подмножество данных о береговой линии из набора данных Глобальной последовательной иерархической географии с высоким разрешением (GSHHG)
Управляйте функциями многоугольника, чтобы добавить озера и другие внутренние водные тела как внутренние звонки многоугольника ("отверстия")
Сохраните модифицированный набор данных в файл форм для будущего использования в MATLAB®, или для экспорта в географическую информационную систему
Глобальная Последовательная Иерархическая География С высоким разрешением (GSHHG; раньше Глобальные Последовательные Иерархические Береговые линии С высоким разрешением или GSHHS), набор данных, Паулем Весселом и Уолтером Х. Ф. Смитом, обеспечивает непротиворечивое множество иерархически расположенных закрытых многоугольников. Они могут использоваться, чтобы создать основные карты, или в приложениях или исследованиях, которые включают операции как географические поисковые запросы или статистические свойства береговых линий.
Этот пример создает несколько временных файлов и использует переменную workingFolder обозначить их местоположение. Значение, используемое здесь, определяется выходом tempdir команда, но вы могли легко настроить это.
workingFolder = tempdir;
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');
Выберите данные для определенного четырехугольника долготы широты и импортируйте их как массив "геоstruct" Mapping Toolbox:
latlim = [-60 15]; lonlim = [-90 -30]; S = gshhs(filename, latlim, lonlim);
Если вы закончили извлекать данные, можно удалить распакованный файл GSHHS и индексный файл.
delete(filename) delete(indexfile)
Исследуйте первый элемент геомассива структур 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 включает четыре уровня береговых линий:
Уровень 1 - "Земля"
Уровень 2 - "Озеро"
Уровень 3 - "Остров в озере"
Уровень 4 - "Водоем на острове в озере"
Проверяйте, чтобы видеть, который выравнивает данные, которые вы импортировали, включает. Level поле содержит числовые числа уровня.
levels = [S.Level]; unique(levels)
ans = 1×3
1 2 3
LevelString поле обеспечивает их интерпретацию. Например,
S(104).LevelString
ans = 'lake'
показывает, что показывают 104, озеро (функция Уровня 2).
В этом примере, должном или к низкому разрешению или к пространственному подмножеству, не присутствуют никакие функции Уровня 4.
Этот пример управляет лучшими двумя уровнями иерархии 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
Задайте анонимную функцию предиката, чтобы обнаружить пересечения ограничительной рамки (возвращающий 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
С одним вызовом shapewrite, можно создать трио файлов,
gshhs_c_SouthAmerica.shp gshhs_c_SouthAmerica.shx gshhs_c_SouthAmerica.dbf
в вашей рабочей папке.
shapepath = fullfile(workingFolder,'gshhs_c_SouthAmerica.shp');
shapewrite(L1, shapepath)Подтверждать результаты 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