Как вычислить градиент матрицы для рисования векторного поля в R?

Question

Я огляделся, но все, что я мог найти, это найти матрицу с использованием diff(d), где d - это матрица. Это не дает мне векторов, просто кучу скаляров. Я не совсем уверен, что с ними делать.

Я хотел бы найти способ вычисления градиента в нескольких точках по всей поверхности, который представлен матрицей. Этот градиент может отображаться как векторное поле. Есть question here о создании векторных полей в R, но я не знаю, как вычислить градиент.

Редактировать: Я попытаюсь подробно остановиться на том, что я ищу. Скажем, у меня есть матрица, как это:

 X0 X1.5 X3.1 X4.3 X5.9 X7.3 X8.6 X9.8 X11 X12.3 X13.6 X14.9 X16.4 X17.9 X20 
[1,] 0 1.4 3.0 4.5 6.0 7.3 8.6 9.7 10.9 12.2 13.4 14.9 16.4 18.1 20 
[2,] 0 1.6 3.2 4.9 6.4 7.6 8.7 9.6 10.6 11.8 13.2 14.7 16.4 18.1 20 
[3,] 0 1.7 3.5 5.2 7.0 8.3 9.0 9.4 9.9 11.1 12.7 14.6 16.3 18.2 20 
[4,] 0 1.8 3.7 5.8 8.0 9.3 9.3 9.3 9.4 10.2 12.1 14.1 16.2 18.0 20 
[5,] 0 1.7 3.9 6.0 8.8 9.3 9.3 9.4 9.6 9.9 11.8 14.0 16.2 18.1 20 
[6,] 0 1.8 3.8 5.7 8.1 9.3 9.3 9.4 9.6 10.1 12.3 14.4 16.3 18.0 20 
[7,] 0 1.6 3.5 5.2 7.0 8.4 9.1 9.5 10.1 11.3 13.0 14.6 16.4 18.2 20 
[8,] 0 1.5 3.2 4.9 6.4 7.7 8.7 9.7 10.7 11.9 13.3 14.9 16.5 18.3 20 
[9,] 0 1.5 3.1 4.6 6.0 7.4 8.6 9.7 10.9 12.1 13.5 15.1 16.6 18.3 20 
[10,] 0 1.5 3.0 4.6 6.0 7.3 8.5 9.7 10.9 12.4 13.6 13.1 16.6 18.2 20 

Это выглядит примерно так, когда вы участка его:

Теперь, что я хочу, это просто: через определенные промежутки времени х и y, я бы хотел найти склон поверхности. Так, например, начиная с x = 0, y = 0, я хотел бы найти наклон в виде вектора, который я могу использовать для построения позже. Затем найдите наклон при x = 0, y = 1 и т. Д. Для всех значений y. Затем найдите все значения y для x = 1 и т. Д.

Цель состоит в том, чтобы иметь связку векторов, которые могут быть построены в векторном поле like this.

Можно ли это сделать в R?

Answer 1

Для начала ознакомьтесь с некоторыми материалами.

m <- matrix(1:9,nrow=3) 

Вы должны решить, нужно ли заполнять НС или 0 в начале или в конце, или повторить первую или последнюю ценность в diff(x), или ...

bdiff <- function(x) c(NA,diff(x)) 

Градиенты в х (строка) направление:

t(apply(m,1,bdiff)) 
##  [,1] [,2] [,3] 
## [1,] NA 3 3 
## [2,] NA 3 3 
## [3,] NA 3 3 

В у (колонка) направлении:

apply(m,2,bdiff) 
##  [,1] [,2] [,3] 
## [1,] NA NA NA 
## [2,] 1 1 1 
## [3,] 1 1 1 

Для примера, что-то примерно так, как это работает:

m2 <- matrix(c(
0,1.4,3.0,4.5,6.0,7.3,8.6,9.7,10.9,12.2,13.4,14.9,16.4,18.1,20, 
0,1.6,3.2,4.9,6.4,7.6,8.7,9.6,10.6,11.8,13.2,14.7,16.4,18.1,20, 
0,1.7,3.5,5.2,7.0,8.3,9.0,9.4,9.9,11.1,12.7,14.6,16.3,18.2,20, 
0,1.8,3.7,5.8,8.0,9.3,9.3,9.3,9.4,10.2,12.1,14.1,16.2,18.0,20, 
0,1.7,3.9,6.0,8.8,9.3,9.3,9.4,9.6,9.9,11.8,14.0,16.2,18.1,20, 
0,1.8,3.8,5.7,8.1,9.3,9.3,9.4,9.6,10.1,12.3,14.4,16.3,18.0,20, 
0,1.6,3.5,5.2,7.0,8.4,9.1,9.5,10.1,11.3,13.0,14.6,16.4,18.2,20, 
0,1.5,3.2,4.9,6.4,7.7,8.7,9.7,10.7,11.9,13.3,14.9,16.5,18.3,20, 
0,1.5,3.1,4.6,6.0,7.4,8.6,9.7,10.9,12.1,13.5,15.1,16.6,18.3,20, 
0,1.5,3.0,4.6,6.0,7.3,8.5,9.7,10.9,12.4,13.6,13.1,16.6,18.2,20), 
byrow=TRUE,nrow=10) 

rr <- row(m2) 
cc <- col(m2) 
dx <- t(apply(m2,1,bdiff)) 
dy <- apply(m2,2,bdiff) 
sc <- 0.25 
off <- -0.5 ## I *think* this is right since we NA'd row=col=1 
plot(rr,cc,col="gray",pch=16) 
arrows(rr+off,cc+off,rr+off+sc*dx,cc+off+sc*dy,length=0.05) 

Answer 2

Вот raster комплексный подход. Начнем с той же матрицей, как ответ Бена:

m2 <- matrix(c(
0,1.4,3.0,4.5,6.0,7.3,8.6,9.7,10.9,12.2,13.4,14.9,16.4,18.1,20, 
0,1.6,3.2,4.9,6.4,7.6,8.7,9.6,10.6,11.8,13.2,14.7,16.4,18.1,20, 
0,1.7,3.5,5.2,7.0,8.3,9.0,9.4,9.9,11.1,12.7,14.6,16.3,18.2,20, 
0,1.8,3.7,5.8,8.0,9.3,9.3,9.3,9.4,10.2,12.1,14.1,16.2,18.0,20, 
0,1.7,3.9,6.0,8.8,9.3,9.3,9.4,9.6,9.9,11.8,14.0,16.2,18.1,20, 
0,1.8,3.8,5.7,8.1,9.3,9.3,9.4,9.6,10.1,12.3,14.4,16.3,18.0,20, 
0,1.6,3.5,5.2,7.0,8.4,9.1,9.5,10.1,11.3,13.0,14.6,16.4,18.2,20, 
0,1.5,3.2,4.9,6.4,7.7,8.7,9.7,10.7,11.9,13.3,14.9,16.5,18.3,20, 
0,1.5,3.1,4.6,6.0,7.4,8.6,9.7,10.9,12.1,13.5,15.1,16.6,18.3,20, 
0,1.5,3.0,4.6,6.0,7.3,8.5,9.7,10.9,12.4,13.6,13.1,16.6,18.2,20), 
byrow=TRUE,nrow=10) 

Преобразовать в растровый (обратите внимание на транспонирование и вообще бесполезный происходит потому, что матрицы начинается с верхнего левого, но координирует работу снизу справа):

require(raster) 
require(rasterVis) 
r=raster(t(m2[,ncol(m2):1]), xmn=0.5,xmx=nrow(m2)+.5, ymn=0.5,ymx=ncol(m2)+0.5) 

Теперь, когда Бен намекнул вам, вы должны дать ему систему координат do. На данный момент он имеет только числовые координаты столбцов и столбцов от 1 до 10 и от 1 до 15. Если это была карта реального мира, то raster должен знать, является ли это длиною в длину или метрами или футами, и будет ли X координаты и координаты Y находятся в одном масштабе.Это важно, даже для данных, которые is not сопоставлены с реальным миром, как я подозреваю, ваши данные.

Градиент не имеет смысла, если ваши координаты X и Y находятся не в тех же единицах. Если X - сопротивление в омах, а Y - ток в амперах, а Z - ваш измеряемый потенциал в вольтах, то какой уклон? Ну, это может быть 2 В на Ом на оси X и -3 В на усилитель в направлении Y. Так что вообще? Вы не можете сказать, потому что вы не можете комбинировать омы и усилители, чтобы получить направление.

Итак, я исхожу из того, что независимо от того, какие единицы X и Y находятся в вашем примере, они являются одинаковыми единицами (возможно, они являются омами на резисторе A и ом на резисторе B), и они идут от 1 до 10 и от 1 до 15.

Теперь я думаю, что есть код проецирования, который просто говорит: «Это координаты x и y без реального географического значения», но я не могу вспомнить, что это такое или найти. Поэтому я просто лежу и использую любую старую систему координат, которую я знаю, является обычной декартовой сеткой. В этом случае - Национальная сеть ГБ. Если вы пытались построить этот растр на карте было бы крошечная площадь у берегов юго-западной Англии, потому что там сетка происхождение, а ваши данные 10м на 15м в этой системе:

projection(r)=CRS("+init=epsg:27700") 

Давайте построить его, чтобы убедиться, что мы не испортили еще:

persp(r,theta=-50,phi=20, shade=0.23,col="red") 

отмечают, что X и координаты Y указывают в том же направлении, что и ваш образец сюжета, так что я знаю, я До сих пор все было в порядке.

Теперь я могу просто сделать levelplot от rasterVis, но мне нужно сделать небольшое масштабирование. Это связано с тем, что градиент на реальной карте вычисляется с высот и расстояний, которые имеют одинаковые единицы (например, метров или футов), но ваши данные являются просто цифрами. Следовательно, градиенты на самом деле довольно малы в натуральной целочисленной системе координат. Итак:

vectorplot(r, scaleSlope=.1) 

дает вам:

Обратите внимание на наклон, как правило, сверху вниз, потому что это путь ваш X и оси Y в вашем примере сюжета (и, следовательно, в моем растр). Также обратите внимание, что ячейки являются квадратными, потому что мы сохраняем соотношение сторон данных (потому что мы обрабатываем координаты X и Y как равные по мере). Ответ Бена показывает общий поток L-R, что означает, что его координаты X и Y не в обычном порядке.

Кроме того, градиент ознакомительного алгоритм vectorplot делает некоторую степень сглаживания, поэтому маленький разрыв в правом верхнем угол не выглядит, как крайним, как в алгоритме разностного Бен:

но вы должны решить, действительно ли вы хотите построить сглаженный градиент или конечные различия ...