Руководство пользователя Upslope Area Toolbox
Содержимое
Что такое "область Upslope"?
Предположите, что вы стоите на стороне выступа где-нибудь во время дождя. Часть воды, которая падает в гору от вашего положения, будет течь непосредственно к и затем мимо вашей обуви. Часть воды, тем не менее, будет течь под гору в различном направлении, далеко от вас. Область земли выше вас, которая высушивает непосредственно через, где вы стоите, называется upslope областью вашего положения.
Если бы вы стояли в самом верху выступа, upslope область, то там был бы 0; никакие потоки воды вам отовсюду еще. С другой стороны, если бы вы стояли в самой глубокой точке в кратере с высокими оправами полностью вокруг, upslope областью была бы целая область кратера.
Областью Upslope является важное измерение гидрологии, используемое, чтобы изучить водные сети дренажа, движение отложений и загрязнителей, эрозии, полных побед.
О цифровых моделях вертикального изменения
Цифровая модель вертикального изменения (или DEM) является компьютерным представлением поверхностной топографии. Растровый DEM является прямолинейной сеткой значений, каждое из которых представляет высоту поверхности в соответствующем местоположении сетки.
Высококачественные, демократы с высоким разрешением теперь широко доступны и используют для большого разнообразия анализа ландшафта. Данные DEM для Соединенных Штатов могут быть получены через американскую Геологическую службу (USGS) и его провайдеров данных.
Файл milford_ma_dem.mat содержит пример DEM, покрывающий фрагмент Массачусетса в Соединенных Штатах. Можно загрузить этот MAT-файл и отобразить матрицу вертикального изменения Z можно следующим образом:
load milford_ma_dem imshow(Z, []) % imshow is in the Image Processing Toolbox
Используйте surf и другие функции графики MATLAB, чтобы отобразить небольшую часть DEM как поверхность.
Zsub = Z(220:250, 170:215); surf(Zsub) shading interp colormap(gray) set(gca, 'YDir', 'reverse') view(-15, 40) axis off
Вы видите яркий выступ на левой стороне и что похоже на нее, может быть водоем (плоская, темная область) в середине. Но воздействие высоты значительно преувеличено здесь. Если вы смотрите на переменную description в MAT-файле, вы видите, что x-и y-разрешение пикселей в данных DEM составляют 30 метров. Можно получить лучшее представление об истинном внешнем виде поверхности путем установки DataAspectRatio соответственно.
set(gca, 'DataAspectRatio', [1 1 30])
view(-10, 10)
Моделирование поверхностных потоков
При анализе потока воды с помощью DEM существенный шаг в анализе должен определить направление потока в каждой точке в сетке DEM. Рассмотрите, например, 3х3 матрицу значений высоты ниже:
E = [10 10.5 11; 10 9 8.9; 10.3 8.5 8.4]
E = 10.0000 10.5000 11.0000 10.0000 9.0000 8.9000 10.3000 8.5000 8.4000
Центральная точка имеет высоту 9. Это - восточный сосед, имеет ту же высоту. Это имеет двух наклонных соседей на юг и юго-восток. Как должен определить направление потока воды для этой точки?
Upslope Area Toolbox обеспечивает функции, которые вычисляют направление потока воды с помощью метода D-бесконечности, описанного в Tarboton, "Новый метод для определения направлений потока и upslope областей в сетке цифровые модели вертикального изменения", Исследование Водных ресурсов, издание 33, № 2, страницы 309-319, февраль 1997.
Функциональный pixelFlow возвращает направление потока для данной точки в DEM. Направление возвращено как угол (в радианах) измеренный против часовой стрелки от указывающей восток горизонтальной оси.
center_point_flow_in_degrees = pixelFlow(E, 2, 2) * (180/pi)
center_point_flow_in_degrees = 281.3099
Таким образом, поток от центральной точки является приблизительно 281 градусом или юго-юго-востоком. Можно использовать demFlow, чтобы вычислить направление потока для всех точек в DEM.
R = demFlow(E)
R =
5.4978 4.7790 4.7124
5.8195 4.9098 4.7124
0 0 NaN
Значение NaN в правом нижнем углу указывает, что местоположение не имеет никаких наклонных соседей, таким образом, нет никакого наклонного потока воды оттуда.
Вычисление и понимание матрицы потока
Другой важный шаг в гидрологическом анализе должен ответить на этот вопрос для каждой точки в DEM: насколько потоки воды в ту точку от каждой из ее соседних точек? Матрица потока, вычисленная функциональным flowMatrix, отвечает на этот вопрос для всех точек в DEM.
T = flowMatrix(E, R)
T = (1,1) 1.0000 (5,1) -1.0000 (2,2) 1.0000 (5,2) -0.4097 (6,2) -0.5903 (3,3) 1.0000 (6,3) -1.0000 (4,4) 1.0000 (5,4) -0.9152 (8,4) -0.0848 (5,5) 1.0000 (6,5) -0.7487 (9,5) -0.2513 (6,6) 1.0000 (9,6) -1.0000 (7,7) 1.0000 (8,7) -1.0000 (8,8) 1.0000 (9,8) -1.0000 (9,9) 1.0000
Несмотря на то, что матрица потока разреженна, в этом небольшом примере легче посмотреть на матричные значения потока, если вы преобразуете его в полный.
T = full(T)
T =
1.0000 0 0 0 0 0 0 0 0
0 1.0000 0 0 0 0 0 0 0
0 0 1.0000 0 0 0 0 0 0
0 0 0 1.0000 0 0 0 0 0
-1.0000 -0.4097 0 -0.9152 1.0000 0 0 0 0
0 -0.5903 -1.0000 0 -0.7487 1.0000 0 0 0
0 0 0 0 0 0 1.0000 0 0
0 0 0 -0.0848 0 0 -1.0000 1.0000 0
0 0 0 0 -0.2513 -1.0000 0 -1.0000 1.0000
Каждая из девяти строк и столбцов T соответствует одной из девяти точек в 3х3 DEM с точками в DEM, пронумерованном по столбцам. Например, четвертая строка и четвертый столбец соответствуют точке DEM E(1,2). Точно так же девятая строка и девятый столбец соответствуют точке DEM E(3,3).
1's по диагонали T представляют идею, что равный единичный объем воды добавляется к поверхности, по-видимому, от дождя, в каждой точке в DEM.
Посмотрите на столбцы T, чтобы видеть, где дождевая вода вытекает в. Вот второй столбец T.
T(:, 2)
ans =
0
1.0000
0
0
-0.4097
-0.5903
0
0
0
Значения в этом столбце указывают, что приблизительно 41% воды, текущей в точку № 2 DEM, течет вниз в точку № 5 DEM, потому что T(5,2) равняется-0.4907. Приблизительно 59% воды, текущей в точку № 2 DEM, текут вниз в точку № 6 DEM, потому что T(6,2) равняется-0.5903.
Посмотрите на строки T, чтобы видеть, где данная точка DEM получает свою воду из. Например, вот девятая строка T:
T(9, :)
ans =
0 0 0 0 -0.2513 -1.0000 0 -1.0000 1.0000
Эти значения указывают, что точка № 9 DEM получает всю воду, текущую в точку № 8 DEM (T(9,8) равняется-1.0), вся вода, текущая в точку № 6 DEM (T(9,6) равняется-1.0), и приблизительно 25% воды, текущей в точку № 5 DEM (T(9,5) равняется-0.2513).
Вычисление область Upslope
Функциональный upslopeArea вычисляет upslope область для каждой точки в DEM путем решения разреженной линейной системы уравнений на основе матрицы потока. Например:
U = upslopeArea(E, T)
U =
1.0000 1.0000 1.0000
1.0000 3.3249 2.0848
1.0000 5.0796 9.0000
Обратите внимание на то, что upslope область точки включает себя в этот расчет. upslope область 1,0 указывает, что единственная вода, текущая в то местоположение, является модульной суммой, принятой от ливня. В нашем небольшом 3х3 примере DEM вода, падающая на все точки в конечном счете, течет под гору в (3,3) местоположение так, чтобы upslope область в той точке равнялась 9.
Теперь давайте решим настоящую задачу. А именно, давайте вычислим upslope область для данных в milford_ma_dem.
load milford_ma_dem
R = demFlow(Z);
T = flowMatrix(Z, R);
U = upslopeArea(Z, T);
imshow(U, [])
Трудно видеть много детали. Другой метод визуализации, который можно попробовать, должен отобразить логарифм upslope области. Этот метод показывает намного больше детали.
imshow(log(U), [])
Можно также использовать функцию visMap, чтобы наложить upslope область (заштрихованный в зеленом) по исходным данным DEM.
visMap(U, Z)
Вот является увеличивший масштаб представлением
axis([165 230 160 290])
Влияйте на карты
Матрица потока может использоваться, чтобы создать другие линейные системы, решения которых дают полезную информацию. Например, можно задать этот вопрос: Для данной точки P DEM, каков полный набор наклонных точек DEM, которые получают воду от P? Карта влияния, матрица, вычисленная influenceMap, отвечает на этот вопрос.
В этом примере вы вычислите карту влияния для milford_ma точки DEM (235, 185) и затем отобразите его с помощью visMap.
I = influenceMap(Z, T, 235, 185);
visMap(I, Z, 235, 185)
% Zoom in
axis([165 230 220 290])
Вы видите, что вода, запускающаяся наверху выступа (синяя точка), течет на восток в водоем и затем через южный конец водоема в локальную переменную минимумы (приемник). (См. раздел "Sinks" в "Специальных Факторах Данных" ниже.)
Карты зависимости
Карта зависимости, матрица, вычисленная dependenceMap, является другим количеством, вычисленным из матрицы потока. Это показывает, что полный набор идущего в гору DEM указывает что дренаж через данное местоположение DEM. Точно так же, как карта влияния карта зависимости может визуализироваться с помощью visMap. Следующий пример показывает, как вычислить и визуализировать карту зависимости для местоположения DEM (270, 189).
D = dependenceMap(Z, T, 270, 189);
visMap(D, Z, 270, 189)
% Zoom in
axis([65 325 70 350])
Существует синяя точка, указывающая на точку, что зеленая зона высушивает через, но в этой шкале это слишком мало, чтобы видеть, таким образом увеличьте масштаб далее:
axis([155 225 235 305])
Специальные факторы данных
Приемники
Для примера карты влияния выше, вспомните, что вода текла вниз от верхней части выступа в водоем, и южный конец водоема, где это, казалось, просто остановилось. Поэтому данные DEM имели локальный минимум там. Вот значения данных DEM, сразу окружающие ту точку.
Z(268:272, 187:191)
ans =
99 98 99 103 104
98 99 99 102 102
100 99 98 99 100
103 102 100 99 99
107 106 103 101 99
Вы видите, что высота в срединной точке равняется 98, который ниже, чем все точки DEM, окружающие его. Этот вид локального минимума называется приемником. Для многих видов топологических исследований желательно устранить все приемники, которые не расположены в ребре DEM. Можно использовать функциональный fillSinks, чтобы устранить эти внутренние приемники.
Давайте повторим пример карты влияния на исходных milford_ma данных DEM, включая это время шаг предварительной обработки, чтобы устранить внутренние приемники.
Zp = fillSinks(Z);
Rp = demFlow(Zp);
Tp = flowMatrix(Zp, Rp);
Ip = influenceMap(Zp, Tp, 235, 185);
visMap(Ip, Zp, 235, 185)
% Zoom in
axis([165 230 220 290])
Теперь вы видите, что вода продолжает течь. Если вы уменьшаете масштаб далее, вы видите, что вода продолжает течь, пока она не достигает минимума, группы равной высоты точек DEM, которая соединяется с ребром набора данных.
axis([100 320 200 473])
Плато
Плато является связанной группой точек DEM с равной высотой. Например, точки DEM в водоеме, показанном здесь, у всех есть значение 106.
imshow(Z, []) axis([170 215 220 250])
Не всегда ясно, как вычислить поток воды через плато. Функциональный demFlow использует arrowing метод, описанный во Ф. Мейере, "Skeletons и Perceptual Graphs", Обработка сигналов 16 (1989) 335-363. (См. Приложение A.2, Arrowing, на страницах 360-361.) Метод работает обоснованно хорошо во многих случаях, но вы видите некоторые артефакты плато в upslope визуализации области для milford_ma DEM.
imshow(log(U), [])
% Zoom in
axis([110 210 20 150])
upslope тулбокс области обеспечивает функциональный postprocessPlateaus, который заменяет upslope значения области для набора точек плато со средним значением upslope область для целого плато. Вот то, как это работало бы на milford_ma данные.
Um = postprocessPlateaus(U, Z);
imshow(log(Um), [])
% Zoom in
axis([110 210 20 150])
Пропавшие без вести Данных, Много наборов данных DEM имеют законные значения только в определенном водоразделе, который обычно является неправильной областью. Массивы DEM обычно имеют "значение заливки" за пределами области водораздела, чтобы указать на недопустимые данные. Обычно значение заливки является чем-то распознаваемым как-999, и обычно значение заливки дано в описании набора данных.
Чтобы обработать такие наборы данных с помощью Upslope Area Toolbox, замените значения заливки на NaN с помощью кода как это:
Z(Z == fill_value) = NaN;
borderNans функции тулбокса идентифицирует NaN-ценные точки DEM, которые соединяются с "внешними" ребрами набора данных DEM. milford_ma данные DEM содержат, ограничивают NaNs, который можно отобразить можно следующим образом:
imshow(borderNans(Z)) % Turn the axis box on so you can see the extent of the white pixels, which are % the border NaNs. axis on
Получение данных DEM
Данные DEM для Соединенных Штатов могут быть получены из американской Геологической службы (USGS). Например, данные DEM и другие наборы данных могут быть получены из Национальной Карты Бесшовный Сервер. Можно загрузить эти данные в формате BIL, который может быть считан с помощью функции MATLAB multibandread. Существует несколько веб-сайтов, которые предлагают информацию и примеры при получении данных из Бесшовного Сервера, включая этого в Йельском университете.
Для определения местоположения данных DEM, покрывающих другие области, вы можете попробовать список "Свободной Цифровой модели вертикального изменения (DEM) и Свободных ссылок на загрузку Спутниковых снимков" в terrainmap.com.
Если у вас есть последняя версия (R2009b или позже) Mapping Toolbox, можно использовать функции картографического Веб-сервиса (WMS), чтобы получить данные DEM. Страница с описанием wmsread имеет пример, показывающий, как получить данные DEM непосредственно из сервера JPL WMS.