Карта эллипсов дисперсии по мат. ожиданию, осям и углу

[nickleb]

Здравствуйте. Кто-нибудь имел опыт в R рисования карты эллипсов дисперсии по мат. ожиданию, осям и углу между его главной осью и мат. ожиданием? Примерно как:
www.pmel.noaa.gov/maillists/tmap/ferret ... 00315.html

[Denis Rykov]

403

[nickleb]

403

Denis, исправил ссылку...

[Максим Дубинин]

исправил ссылку

нет

[nickleb]

исправил ссылку

нет

извините... вот так сейчас: и здесь, и в первом сообщении темы:
www.pmel.noaa.gov/maillists/tmap/ferret ... 00315.html

[gamm]

Код: Выделить всё

my.ellipse<-function(x0,y0,a,b,ell.angle,ell.segm=48,print=TRUE) {
  # horisonl ellipse
  xy<-matrix(ncol=2,nrow=ell.segm+1)
  for(i in 1:ell.segm) {
    ang<-2*pi/ell.segm
    x<-cos((i-1)*ang)*a
    y<-sin((i-1)*ang)*b
    xy[i,]<-c(x,y)
  }
  xy[ell.segm+1,]<-xy[1,]
  
  # turn ellipse
  ell.sin<-sin(ell.angle/180*pi)
  ell.cos<-cos(ell.angle/180*pi)
  mtr<-matrix(c(ell.cos,ell.sin,-ell.sin,ell.cos),ncol=2,nrow=2)
  xy1<-t((mtr %*% t(xy)) + c(x0,y0))
  if(print) points(xy1,type="l",lwd=3)
  invisible(xy1)
}

pos<-expand.grid(1:10*10,1:5*10)
plot(0,0,xlim=c(0,105),ylim=c(0,55),asp=1,type="n")
for(i in 1:nrow(pos)) {
  x0<-pos[i,1]        # center position 
  y0<-pos[i,2]        # center position 
  a<-runif(1,0.1,5)   # "main" semiaxis
  b<-runif(1,0.1,5)   # "secondary" semiaxis
  ang<-runif(1,0,360) # angle
  my.ellipse(x0,y0,a,b,ang,ell.segm=48,print=TRUE)
}

главная полуось горизонтальна, поворот - против часовой от горизонтали.

[nickleb]

Код: Выделить всё

my.ellipse<-function(x0,y0,a,b,ell.angle,ell.segm=48,print=TRUE) {
  # horisonl ellipse
  xy<-matrix(ncol=2,nrow=ell.segm+1)
  for(i in 1:ell.segm) {
    ang<-2*pi/ell.segm
    x<-cos((i-1)*ang)*a
    y<-sin((i-1)*ang)*b
    xy[i,]<-c(x,y)
  }
  xy[ell.segm+1,]<-xy[1,]
  
  # turn ellipse
  ell.sin<-sin(ell.angle/180*pi)
  ell.cos<-cos(ell.angle/180*pi)
  mtr<-matrix(c(ell.cos,ell.sin,-ell.sin,ell.cos),ncol=2,nrow=2)
  xy1<-t((mtr %*% t(xy)) + c(x0,y0))
  if(print) points(xy1,type="l",lwd=3)
  invisible(xy1)
}

pos<-expand.grid(1:10*10,1:5*10)
plot(0,0,xlim=c(0,105),ylim=c(0,55),asp=1,type="n")
for(i in 1:nrow(pos)) {
  x0<-pos[i,1]        # center position 
  y0<-pos[i,2]        # center position 
  a<-runif(1,0.1,5)   # "main" semiaxis
  b<-runif(1,0.1,5)   # "secondary" semiaxis
  ang<-runif(1,0,360) # angle
  my.ellipse(x0,y0,a,b,ang,ell.segm=48,print=TRUE)
}

главная полуось горизонтальна, поворот - против часовой от горизонтали.

Ваш пример сработал, gamm. Спасибо! Буду пробовать адаптировать для себя.