研究のページ/R/Rでインタラクティブなグラフを作る
をテンプレートにして作成
[
トップ
] [
新規
| |
検索
|
最終更新
|
ヘルプ
|
ログイン
]
開始行:
#topicpath
* Rでインタラクティブなグラフを作る [#m5e4133f]
** 参考URL [#md84b646]
- 英語
-- [[Plotly R>https://plotly.com/r/]]
-- [[Interactive web-based data visualization with R, plo...
-- [[Scatter plot in R using Plotly>http://rstudio-pubs-s...
-- [[htmlwidgets for R>http://www.htmlwidgets.org/]]
- 日本語
-- [[Kazutan.R / plotly入門>https://kazutan.github.io/kaz...
-- [[plotlyと組み合わせてShinyアプリケーションを作成する>...
** 実際に作った例 [#q3d2a7e7]
*** 有限振幅振り子のハミルトニアンの3D描画 [#k8e0d1e6]
+ コード
#divregion
##########################################################
## plotly で 有限振幅振り子の Hamiltonian
##########################################################
## constants
m <- 1
g <- 10
l <- 0.5
## make plot data
library(plotly)
axx <- list( title = "position q")
axy <- list( title = "momentum p")
axz <- list( title = "Hamiltonian H(q,p)")
x <- seq( -3*pi, 3*pi, pi/30)
y <- seq( -4.7, 4.7, 0.1)
z <- matrix( rep( 0, length(x) * length(y) ), length(y))
for( i in 1:length(y))
for( j in 1:length(x)){
z[i,j] <- y[i]*y[i]/2/m + m*g*l*(1-cos(x[j]))
}
fig <- plot_ly(x=~x, y=~y, z = ~z)
fig <- fig %>% layout(scene = list(xaxis=axx,yaxis=axy,...
fig <- fig %>% add_surface(
contours = list(
z = list(
show=TRUE,
usecolormap=TRUE,
highlightcolor="#ff0000",
project=list(z=TRUE)
)
)
)
fig
#enddivregion
+ [[出来上がったもの>https://robo.mydns.jp/Lecture/HTML/R...
#iframe(https://robo.mydns.jp/Lecture/HTML/R_Plotly/viewh...
*** 球面調和関数 [#x99f208b]
+ コード~
#divregion
########################################################...
## plot_ly の mesh3d で球面上の分布を描くためのメッシュ...
## cf. https://plotly.com/r/reference/#mesh3d
## cf. https://plotly.com/r/3d-mesh/
########################################################...
## 返り値
## x, y, z グリッドの座標
## vi, vj, vk 三角形の頂点を、上述のx,y,z の何番目に...
## ポイントの番号はゼロから始まることに注意。
########################################################...
library(plotly)
spherical3dmesh <- function(nth=21, nph=40){
## θとφの決定
theta <- seq(pi/(nth-1), pi-pi/(nth-1), pi/(nth-1) );
phi <- seq(0, 2*pi-pi/nth, 2*pi/nph );
## x, y, z の決定
x <- c(0)
y <- c(0)
z <- c(1)
for( i in 1:(nth-2) ){
for( j in 1:nph ){
x <- c(x, sin(theta[i]) * cos(phi[j]))
y <- c(y, sin(theta[i]) * sin(phi[j]))
z <- c(z, cos(theta[i]) )
}
}
x <- c(x, 0)
y <- c(y, 0)
z <- c(z,-1)
## vi,vj,vz (三角形の頂点)の決定
### 北極を囲む三角形
vi <- c()
vj <- c()
vk <- c()
for( i in 1:(nph-1) ){
vi <- c(vi,0)
vj <- c(vj,i)
vk <- c(vk,1+i)
}
vi <- c(vi, 0)
vj <- c(vj, rev(vk)[1])
vk <- c(vk, vj[1])
### 側面の三角形
for( i in 1:(nth-3) ){
startn <- rev(vk)[1]
starts <- rev(vj)[1]+1
for( j in 1:(nph-1) ){
vi <- c(vi, startn+j-1)
vj <- c(vj, starts+j-1)
vk <- c(vk, startn+j)
vi <- c(vi, startn+j)
vj <- c(vj, starts+j-1)
vk <- c(vk, starts+j)
}
endn <- rev(vk)[2]
ends <- rev(vk)[1]
vi <- c(vi, endn)
vj <- c(vj, ends)
vk <- c(vk, startn)
vi <- c(vi, startn)
vj <- c(vj, ends)
vk <- c(vk, starts)
}
### 南極を囲む三角形
startn <- rev(vk)[1]
starts <- length(x)-1
for( j in 1:(nph-1) ){
vi <- c(vi,startn+j)
vj <- c(vj,startn+j-1)
vk <- c(vk,starts)
}
vj <- c(vj,rev(vi)[1])
vi <- c(vi,startn)
vk <- c(vk,starts)
return( list(x=x, y=y, z=z, vi=vi, vj=vj, vk=vk))
}
## ルジャンドル多項式・ルジャンドル陪多項式
## もっと、かっこよく書くのは将来の課題。
Legendre <- function( m, l, z ){
if( m>l ) return( 0 * z )
if( m==0 ){
if( l == 0 ) return( 0 * z + 1/2/sqrt(pi))
else if( l == 1 ) return( 1/2 * sqrt(3/pi) * z )
else if( l == 2 ) return( 1/4 * sqrt(5/pi) * (3 ...
else if( l == 3 ) return( 1/4 * sqrt(7/pi) * (5 ...
}else if( m== 1){
if( l == 1 ) return(-1/2 * sqrt( 3/2/pi) *...
else if( l == 2 ) return(-1/2 * sqrt( 15/2/pi) *...
else if( l == 3 ) return(-1/8 * sqrt( 21/ pi) *...
}else if( m== 2){
if( l == 2 ) return( 1/4 * sqrt( 15/2/pi) *...
else if( l == 3 ) return( 1/4 * sqrt(105/2/pi) *...
}else if( m== 3){
return(-1/8 * sqrt( 35/ pi) *...
}
}
########################################################...
## 本体部分
points <- spherical3dmesh()
m <- 2
l <- 3
fig <- plot_ly(type='mesh3d',
x=points$x,y=points$y,z=points$z,i=points...
## intensity=3*points$x
intensity= cos( m * atan2(points$y, point...
)
htmlwidgets::saveWidget(as_widget(fig), paste("M",m,"L",...
#enddivregion
+ [[出来上がったもの>https://robo.mydns.jp/Lecture/HTML/R...
#iframe(https://robo.mydns.jp/Lecture/HTML/R_Plotly/SPHER...
*** NCEP/NCAR のデータの球面上へのプロット [#i82ca370]
+ コード
-- 本体
#divregion
library(plotly)
library(ncdf4)
library(globe)
library(lubridate)
source("00_Functions.R")
###############################################
# DATA : vwnd
ncin <- nc_open("vwnd.2022.nc")
dname <- "vwnd"
vwnd <- ncvar_get(ncin,dname)
str(vwnd)
lon <- ncvar_get(ncin,"lon")
lat <- ncvar_get(ncin,"lat")
time <- ncvar_get(ncin,"time")
level <- ncvar_get(ncin,"level")
tunits <- ncatt_get(ncin,"time","units")
points <- spherical3dmesh( 73, 144 )
fig <- plot_ly(
x=points$x,y=points$y,z=points$z,i=points...
type='mesh3d',
intensity= c( 0, (vwnd[,,10,940])[145:103...
cmin=-45,
cmax= 45,
colors = c("#8888FF", "#FFFFFF", "#FF8888")
)
fig <- sphericalmap(fig) %>%
layout(title=paste( ymd_hms("1800-01-01 00:00:00") +...
htmlwidgets::saveWidget(as_widget(fig), "vwnd_200hPa_940...
#enddivregion
-- 関数
#divregion
########################################################...
## plot_ly の mesh3d で球面上の分布を描くためのメッシュ...
## cf. https://plotly.com/r/reference/#mesh3d
## cf. https://plotly.com/r/3d-mesh/
########################################################...
## 返り値
## x, y, z グリッドの座標
## vi, vj, vk 三角形の頂点を、上述のx,y,z の何番目に...
## ポイントの番号はゼロから始まることに注意。
########################################################...
spherical3dmesh <- function(nth=21, nph=40){
## θとφの決定
theta <- seq(pi/(nth-1), pi-pi/(nth-1), pi/(nth-1) );
phi <- seq(0, 2*pi-pi/nth, 2*pi/nph );
## x, y, z の決定
x <- c(0)
y <- c(0)
z <- c(1)
for( i in 1:(nth-2) ){
for( j in 1:nph ){
x <- c(x, sin(theta[i]) * cos(phi[j]))
y <- c(y, sin(theta[i]) * sin(phi[j]))
z <- c(z, cos(theta[i]) )
}
}
x <- c(x, 0)
y <- c(y, 0)
z <- c(z,-1)
## vi,vj,vz (三角形の頂点)の決定
### 北極を囲む三角形
vi <- c()
vj <- c()
vk <- c()
for( i in 1:(nph-1) ){
vi <- c(vi,0)
vj <- c(vj,i)
vk <- c(vk,1+i)
}
vi <- c(vi, 0)
vj <- c(vj, rev(vk)[1])
vk <- c(vk, vj[1])
### 側面の三角形
for( i in 1:(nth-3) ){
startn <- rev(vk)[1]
starts <- rev(vj)[1]+1
for( j in 1:(nph-1) ){
vi <- c(vi, startn+j-1)
vj <- c(vj, starts+j-1)
vk <- c(vk, startn+j)
vi <- c(vi, startn+j)
vj <- c(vj, starts+j-1)
vk <- c(vk, starts+j)
}
endn <- rev(vk)[2]
ends <- rev(vk)[1]
vi <- c(vi, endn)
vj <- c(vj, ends)
vk <- c(vk, startn)
vi <- c(vi, startn)
vj <- c(vj, ends)
vk <- c(vk, starts)
}
### 南極を囲む三角形
startn <- rev(vk)[1]
starts <- length(x)-1
for( j in 1:(nph-1) ){
vi <- c(vi,startn+j)
vj <- c(vj,startn+j-1)
vk <- c(vk,starts)
}
vj <- c(vj,rev(vi)[1])
vi <- c(vi,startn)
vk <- c(vk,starts)
return( list(x=x, y=y, z=z, vi=vi, vj=vj, vk=vk))
}
########################################################...
##
## 世界地図の海岸線を得るための関数(globe版)
##
globedata3dxyz <- function(r=1.0){
data("earth")
x <- r * cos( earth$coords[,2] * pi / 180 ) * cos( ...
y <- r * cos( earth$coords[,2] * pi / 180 ) * sin( ...
z <- r * sin( earth$coords[,2] * pi / 180 )
result <- NULL
count <- 1
islands <- 1
while( count < length(x) ){
result[[islands]] <- data.frame(x=x[count:(-1+c...
y=y[count:(-1+c...
z=z[count:(-1+c...
# print("end")
# print(-1+count+earth$runlen[islands+1])
count <- count+earth$runlen[islands+1]
# print(count)
islands <- islands +1
}
return(list(max=islands-1,data=result))
}
sphericalmap <- function(fig){
# mapdata3d <- getworldmap3dxyz()
mapdata3d <- globedata3dxyz()
for( i in 1:mapdata3d$max ){
# for( i in 1:30 ){
fig <- add_trace(fig,
type='scatter3d', mode='lines',
data=as.data.frame(mapdata3d[[...
x=~x, y=~y, z=~z,
i=NULL, j=NULL, k=NULL, intens...
name="",
ids=NULL,
showlegend=FALSE,
line=list(width=1, color='gray...
}
return(fig)
}
#enddivregion
+ [[出来上がったもの>https://robo.mydns.jp/Lecture/HTML/R...
※ 球面に地図と図を描くのは、demo でも見ることができる。
demo("sf-plotly-3D-globe", package = "plotly")
#iframe(https://robo.mydns.jp/Lecture/HTML/R_Plotly/NCEPN...
終了行:
#topicpath
* Rでインタラクティブなグラフを作る [#m5e4133f]
** 参考URL [#md84b646]
- 英語
-- [[Plotly R>https://plotly.com/r/]]
-- [[Interactive web-based data visualization with R, plo...
-- [[Scatter plot in R using Plotly>http://rstudio-pubs-s...
-- [[htmlwidgets for R>http://www.htmlwidgets.org/]]
- 日本語
-- [[Kazutan.R / plotly入門>https://kazutan.github.io/kaz...
-- [[plotlyと組み合わせてShinyアプリケーションを作成する>...
** 実際に作った例 [#q3d2a7e7]
*** 有限振幅振り子のハミルトニアンの3D描画 [#k8e0d1e6]
+ コード
#divregion
##########################################################
## plotly で 有限振幅振り子の Hamiltonian
##########################################################
## constants
m <- 1
g <- 10
l <- 0.5
## make plot data
library(plotly)
axx <- list( title = "position q")
axy <- list( title = "momentum p")
axz <- list( title = "Hamiltonian H(q,p)")
x <- seq( -3*pi, 3*pi, pi/30)
y <- seq( -4.7, 4.7, 0.1)
z <- matrix( rep( 0, length(x) * length(y) ), length(y))
for( i in 1:length(y))
for( j in 1:length(x)){
z[i,j] <- y[i]*y[i]/2/m + m*g*l*(1-cos(x[j]))
}
fig <- plot_ly(x=~x, y=~y, z = ~z)
fig <- fig %>% layout(scene = list(xaxis=axx,yaxis=axy,...
fig <- fig %>% add_surface(
contours = list(
z = list(
show=TRUE,
usecolormap=TRUE,
highlightcolor="#ff0000",
project=list(z=TRUE)
)
)
)
fig
#enddivregion
+ [[出来上がったもの>https://robo.mydns.jp/Lecture/HTML/R...
#iframe(https://robo.mydns.jp/Lecture/HTML/R_Plotly/viewh...
*** 球面調和関数 [#x99f208b]
+ コード~
#divregion
########################################################...
## plot_ly の mesh3d で球面上の分布を描くためのメッシュ...
## cf. https://plotly.com/r/reference/#mesh3d
## cf. https://plotly.com/r/3d-mesh/
########################################################...
## 返り値
## x, y, z グリッドの座標
## vi, vj, vk 三角形の頂点を、上述のx,y,z の何番目に...
## ポイントの番号はゼロから始まることに注意。
########################################################...
library(plotly)
spherical3dmesh <- function(nth=21, nph=40){
## θとφの決定
theta <- seq(pi/(nth-1), pi-pi/(nth-1), pi/(nth-1) );
phi <- seq(0, 2*pi-pi/nth, 2*pi/nph );
## x, y, z の決定
x <- c(0)
y <- c(0)
z <- c(1)
for( i in 1:(nth-2) ){
for( j in 1:nph ){
x <- c(x, sin(theta[i]) * cos(phi[j]))
y <- c(y, sin(theta[i]) * sin(phi[j]))
z <- c(z, cos(theta[i]) )
}
}
x <- c(x, 0)
y <- c(y, 0)
z <- c(z,-1)
## vi,vj,vz (三角形の頂点)の決定
### 北極を囲む三角形
vi <- c()
vj <- c()
vk <- c()
for( i in 1:(nph-1) ){
vi <- c(vi,0)
vj <- c(vj,i)
vk <- c(vk,1+i)
}
vi <- c(vi, 0)
vj <- c(vj, rev(vk)[1])
vk <- c(vk, vj[1])
### 側面の三角形
for( i in 1:(nth-3) ){
startn <- rev(vk)[1]
starts <- rev(vj)[1]+1
for( j in 1:(nph-1) ){
vi <- c(vi, startn+j-1)
vj <- c(vj, starts+j-1)
vk <- c(vk, startn+j)
vi <- c(vi, startn+j)
vj <- c(vj, starts+j-1)
vk <- c(vk, starts+j)
}
endn <- rev(vk)[2]
ends <- rev(vk)[1]
vi <- c(vi, endn)
vj <- c(vj, ends)
vk <- c(vk, startn)
vi <- c(vi, startn)
vj <- c(vj, ends)
vk <- c(vk, starts)
}
### 南極を囲む三角形
startn <- rev(vk)[1]
starts <- length(x)-1
for( j in 1:(nph-1) ){
vi <- c(vi,startn+j)
vj <- c(vj,startn+j-1)
vk <- c(vk,starts)
}
vj <- c(vj,rev(vi)[1])
vi <- c(vi,startn)
vk <- c(vk,starts)
return( list(x=x, y=y, z=z, vi=vi, vj=vj, vk=vk))
}
## ルジャンドル多項式・ルジャンドル陪多項式
## もっと、かっこよく書くのは将来の課題。
Legendre <- function( m, l, z ){
if( m>l ) return( 0 * z )
if( m==0 ){
if( l == 0 ) return( 0 * z + 1/2/sqrt(pi))
else if( l == 1 ) return( 1/2 * sqrt(3/pi) * z )
else if( l == 2 ) return( 1/4 * sqrt(5/pi) * (3 ...
else if( l == 3 ) return( 1/4 * sqrt(7/pi) * (5 ...
}else if( m== 1){
if( l == 1 ) return(-1/2 * sqrt( 3/2/pi) *...
else if( l == 2 ) return(-1/2 * sqrt( 15/2/pi) *...
else if( l == 3 ) return(-1/8 * sqrt( 21/ pi) *...
}else if( m== 2){
if( l == 2 ) return( 1/4 * sqrt( 15/2/pi) *...
else if( l == 3 ) return( 1/4 * sqrt(105/2/pi) *...
}else if( m== 3){
return(-1/8 * sqrt( 35/ pi) *...
}
}
########################################################...
## 本体部分
points <- spherical3dmesh()
m <- 2
l <- 3
fig <- plot_ly(type='mesh3d',
x=points$x,y=points$y,z=points$z,i=points...
## intensity=3*points$x
intensity= cos( m * atan2(points$y, point...
)
htmlwidgets::saveWidget(as_widget(fig), paste("M",m,"L",...
#enddivregion
+ [[出来上がったもの>https://robo.mydns.jp/Lecture/HTML/R...
#iframe(https://robo.mydns.jp/Lecture/HTML/R_Plotly/SPHER...
*** NCEP/NCAR のデータの球面上へのプロット [#i82ca370]
+ コード
-- 本体
#divregion
library(plotly)
library(ncdf4)
library(globe)
library(lubridate)
source("00_Functions.R")
###############################################
# DATA : vwnd
ncin <- nc_open("vwnd.2022.nc")
dname <- "vwnd"
vwnd <- ncvar_get(ncin,dname)
str(vwnd)
lon <- ncvar_get(ncin,"lon")
lat <- ncvar_get(ncin,"lat")
time <- ncvar_get(ncin,"time")
level <- ncvar_get(ncin,"level")
tunits <- ncatt_get(ncin,"time","units")
points <- spherical3dmesh( 73, 144 )
fig <- plot_ly(
x=points$x,y=points$y,z=points$z,i=points...
type='mesh3d',
intensity= c( 0, (vwnd[,,10,940])[145:103...
cmin=-45,
cmax= 45,
colors = c("#8888FF", "#FFFFFF", "#FF8888")
)
fig <- sphericalmap(fig) %>%
layout(title=paste( ymd_hms("1800-01-01 00:00:00") +...
htmlwidgets::saveWidget(as_widget(fig), "vwnd_200hPa_940...
#enddivregion
-- 関数
#divregion
########################################################...
## plot_ly の mesh3d で球面上の分布を描くためのメッシュ...
## cf. https://plotly.com/r/reference/#mesh3d
## cf. https://plotly.com/r/3d-mesh/
########################################################...
## 返り値
## x, y, z グリッドの座標
## vi, vj, vk 三角形の頂点を、上述のx,y,z の何番目に...
## ポイントの番号はゼロから始まることに注意。
########################################################...
spherical3dmesh <- function(nth=21, nph=40){
## θとφの決定
theta <- seq(pi/(nth-1), pi-pi/(nth-1), pi/(nth-1) );
phi <- seq(0, 2*pi-pi/nth, 2*pi/nph );
## x, y, z の決定
x <- c(0)
y <- c(0)
z <- c(1)
for( i in 1:(nth-2) ){
for( j in 1:nph ){
x <- c(x, sin(theta[i]) * cos(phi[j]))
y <- c(y, sin(theta[i]) * sin(phi[j]))
z <- c(z, cos(theta[i]) )
}
}
x <- c(x, 0)
y <- c(y, 0)
z <- c(z,-1)
## vi,vj,vz (三角形の頂点)の決定
### 北極を囲む三角形
vi <- c()
vj <- c()
vk <- c()
for( i in 1:(nph-1) ){
vi <- c(vi,0)
vj <- c(vj,i)
vk <- c(vk,1+i)
}
vi <- c(vi, 0)
vj <- c(vj, rev(vk)[1])
vk <- c(vk, vj[1])
### 側面の三角形
for( i in 1:(nth-3) ){
startn <- rev(vk)[1]
starts <- rev(vj)[1]+1
for( j in 1:(nph-1) ){
vi <- c(vi, startn+j-1)
vj <- c(vj, starts+j-1)
vk <- c(vk, startn+j)
vi <- c(vi, startn+j)
vj <- c(vj, starts+j-1)
vk <- c(vk, starts+j)
}
endn <- rev(vk)[2]
ends <- rev(vk)[1]
vi <- c(vi, endn)
vj <- c(vj, ends)
vk <- c(vk, startn)
vi <- c(vi, startn)
vj <- c(vj, ends)
vk <- c(vk, starts)
}
### 南極を囲む三角形
startn <- rev(vk)[1]
starts <- length(x)-1
for( j in 1:(nph-1) ){
vi <- c(vi,startn+j)
vj <- c(vj,startn+j-1)
vk <- c(vk,starts)
}
vj <- c(vj,rev(vi)[1])
vi <- c(vi,startn)
vk <- c(vk,starts)
return( list(x=x, y=y, z=z, vi=vi, vj=vj, vk=vk))
}
########################################################...
##
## 世界地図の海岸線を得るための関数(globe版)
##
globedata3dxyz <- function(r=1.0){
data("earth")
x <- r * cos( earth$coords[,2] * pi / 180 ) * cos( ...
y <- r * cos( earth$coords[,2] * pi / 180 ) * sin( ...
z <- r * sin( earth$coords[,2] * pi / 180 )
result <- NULL
count <- 1
islands <- 1
while( count < length(x) ){
result[[islands]] <- data.frame(x=x[count:(-1+c...
y=y[count:(-1+c...
z=z[count:(-1+c...
# print("end")
# print(-1+count+earth$runlen[islands+1])
count <- count+earth$runlen[islands+1]
# print(count)
islands <- islands +1
}
return(list(max=islands-1,data=result))
}
sphericalmap <- function(fig){
# mapdata3d <- getworldmap3dxyz()
mapdata3d <- globedata3dxyz()
for( i in 1:mapdata3d$max ){
# for( i in 1:30 ){
fig <- add_trace(fig,
type='scatter3d', mode='lines',
data=as.data.frame(mapdata3d[[...
x=~x, y=~y, z=~z,
i=NULL, j=NULL, k=NULL, intens...
name="",
ids=NULL,
showlegend=FALSE,
line=list(width=1, color='gray...
}
return(fig)
}
#enddivregion
+ [[出来上がったもの>https://robo.mydns.jp/Lecture/HTML/R...
※ 球面に地図と図を描くのは、demo でも見ることができる。
demo("sf-plotly-3D-globe", package = "plotly")
#iframe(https://robo.mydns.jp/Lecture/HTML/R_Plotly/NCEPN...
ページ名: