В этом примере показано, как к:
Извлеките подмножество данных о береговой линии из набора данных Глобальной последовательной иерархической географии с высоким разрешением (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: [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 включает четыре уровня береговых линий:
Уровень 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');
shapewrite(L1, shapepath)
Подтверждать результаты shapewrite
, считайте новый файл форм в геомассив структур southAmerica
:
southAmerica = shaperead(shapepath, 'UseGeoCoords', true)
southAmerica=79×16 struct
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
gshhs
| shaperead
| shapewrite