![[PukiWiki]](image/ROBO192.png "[PukiWiki]")
#author("2022-09-11T19:28:41+09:00","external:moriat","moriat") #author("2022-09-12T11:27:52+09:00","external:moriat","moriat") #topicpath * Rでインタラクティブなグラフを作る [#m5e4133f] ** 参考URL [#md84b646] - 英語 -- [[Plotly R>https://plotly.com/r/]] -- [[Interactive web-based data visualization with R, plotly, and shiny(書籍の電子版)>https://plotly-r.com/]] -- [[Scatter plot in R using Plotly>http://rstudio-pubs-static.s3.amazonaws.com/502499_d17e66bdecfa4e698ea1f8b3fe379395.html]] -- [[htmlwidgets for R>http://www.htmlwidgets.org/]] - 日本語 -- [[Kazutan.R / plotly入門>https://kazutan.github.io/kazutanR/plotly_intro.html]] -- [[plotlyと組み合わせてShinyアプリケーションを作成する>https://www.randpy.tokyo/entry/shiny_31]] ** 実際に作った例 [#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,zaxis=axz)) 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_Plotly/viewhtml3e093b8a7e93/index.html]]~ #iframe(https://robo.mydns.jp/Lecture/HTML/R_Plotly/viewhtml3e093b8a7e93/index.html,style=width:800px;height:480px;)) *** 球面調和関数 [#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 * z^2 - 1) ) else if( l == 3 ) return( 1/4 * sqrt(7/pi) * (5 * z^3 - 3 * z) ) }else if( m== 1){ if( l == 1 ) return(-1/2 * sqrt( 3/2/pi) * sqrt(1-z^2) ) else if( l == 2 ) return(-1/2 * sqrt( 15/2/pi) * sqrt(1-z^2)*z ) else if( l == 3 ) return(-1/8 * sqrt( 21/ pi) * sqrt(1-z^2)*(5 * z^2 - 1) ) }else if( m== 2){ if( l == 2 ) return( 1/4 * sqrt( 15/2/pi) * (1-z^2) ) else if( l == 3 ) return( 1/4 * sqrt(105/2/pi) * (1-z^2)*z ) }else if( m== 3){ return(-1/8 * sqrt( 35/ pi) * sqrt(1-z^2)^3) } } ############################################################################################## ## 本体部分 points <- spherical3dmesh() m <- 2 l <- 3 fig <- plot_ly(type='mesh3d', x=points$x,y=points$y,z=points$z,i=points$vi,j=points$vj,k=points$vk, ## intensity=3*points$x intensity= cos( m * atan2(points$y, points$x)) * Legendre(m,l,points$z) ) htmlwidgets::saveWidget(as_widget(fig), paste("M",m,"L",l,".html", sep="")) #enddivregion + [[出来上がったもの>https://robo.mydns.jp/Lecture/HTML/R_Plotly/SPHERICAL/M2L3.html]]~ #iframe(https://robo.mydns.jp/Lecture/HTML/R_Plotly/SPHERICAL/M2L3.html,style=width:800px;height:480px;)) *** 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$vi,j=points$vj,k=points$vk, type='mesh3d', intensity= c( 0, (vwnd[,,10,940])[145:10368], 0 ), cmin=-45, cmax= 45, colors = c("#8888FF", "#FFFFFF", "#FF8888") ) fig <- sphericalmap(fig) %>% layout(title=paste( ymd_hms("1800-01-01 00:00:00") + hours(time[940]), "Wind V" ) ) htmlwidgets::saveWidget(as_widget(fig), "vwnd_200hPa_940.html") #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( earth$coords[,1] * pi / 180 ) y <- r * cos( earth$coords[,2] * pi / 180 ) * sin( earth$coords[,1] * pi / 180 ) 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+count+earth$runlen[islands+1])], y=y[count:(-1+count+earth$runlen[islands+1])], z=z[count:(-1+count+earth$runlen[islands+1])]) # 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[[2]][[i]]), x=~x, y=~y, z=~z, i=NULL, j=NULL, k=NULL, intensity=NULL, cmin=NULL, cmax=NULL, name="", ids=NULL, showlegend=FALSE, line=list(width=1, color='gray')) } return(fig) } #enddivregion + [[出来上がったもの>https://robo.mydns.jp/Lecture/HTML/R_Plotly/NCEPNCARRe/vwnd_200hPa_940.html]]~ ※ 凡例(legend)がうまくいっていない。~ ※ scatter3d と layout.geo との組み合わせは、demo でも見ることができる。 ※ 球面に地図と図を描くのは、demo でも見ることができる。 demo("sf-plotly-3D-globe", package = "plotly") #iframe(https://robo.mydns.jp/Lecture/HTML/R_Plotly/NCEPNCARRe/vwnd_200hPa_940.html,style=width:800px;height:480px;))
