亡八の覚書

環境分野の化学分析屋(見習い)が、独学でたどる解析の冥府魔道

信頼区間95%が真値を含む様子の図

データの95%信頼区間

場合によっては真値を含まないケースを図示するので

以前、どこかを参考に書いたはずの絵が見つからないので

ggplot2で描いてみた。

もっと良い方法があるのかもしれないけど

ggplot2の練習もかねて、これで一応、良しとする。

 

dataset <- NULL
setmean <- 5
for (i in 1:100){
sampledata <- rnorm(100, setmean, 1)
testresu <- t.test(sampledata)
check <- NULL
if (testresu$conf.int[2] < setmean | testresu$conf.int[1] > setmean){
check <- 1
} else{
check <- 2
}
dataset <- rbind(dataset, c(i, testresu$estimate, testresu$conf.int[1], testresu$conf.int[2], check))
}
colnames(dataset) <- c("No", "ave", "min", "max", "check")
dataset <- as.data.frame(dataset)


ggplot() +
geom_line(aes(x=c(0,101), y =c(setmean, setmean)), size=2, colour="#00FF00") +
geom_point(data = subset(dataset, dataset$check == 1), aes(x=No, y=ave),
colour = "#FF0000") +
geom_errorbar(data = subset(dataset, dataset$check == 1),
aes(x=No, ymin = min, ymax = max), width =.2,
colour = "#FF0000") +
geom_point(data = subset(dataset, dataset$check == 2), aes(x=No, y=ave)) +
geom_errorbar(data = subset(dataset, dataset$check == 2),
aes(x=No, ymin = min, ymax = max), width =.2,
colour = "#000000") +
xlab("") + ylab("") +
theme(plot.title = element_text(size=20),
legend.title = element_text(size=15),
legend.text = element_text(size=15),
axis.text.x=element_blank(),
axis.text.y=element_text(size=15),
axis.title.y = element_text(size=15),
panel.border= element_rect(colour="#000000", fill=NA, size=1.5),
panel.background = element_rect(fill="#FFFFFF"),
panel.grid.major = element_blank())

 

f:id:Q21:20171206141817p:plain

GMTの絵をPSファイル以外に

GMTで書いたものをPSファイル以外にしたいと

思っていたら、ありました。

 

単なる勉強不足でしたなぁ。

 

psconvert map.ps -Tg -E250 -A

 

psconvertを使えばいいとか。

-TgはPNGにするオプション。

-E250は解像度。

-Aは余白を極力削るオプション。

 

まだまだオプションはありそうだけど

大体こんな感じで出来た。

 

覚えておかないと。

R:検索結果がヒットしないのに、NAとか0にはならないとき

Rでのちょっと戸惑ったこと。

 

データからsubset()で検索値をとりだして

その値を使ってさらに足し算したいとき。

検索結果がヒットしなかったら

0を返してほしかったのだけど

そうじゃなくてlogical(0)なんて出てきて

これが煮ても焼いても食えなかった。

 

is.na()をやってもだめで

この数字自体に1をかけたり、1を割ったりして

何かエラーにしようとしても

numeric(0)とか出て、うまく扱えない。

 

で、length()で大きさを見ると、見事に数字の0に落ち着いた。

 

d <- subset(data, lon == basedata$lon[ i ] & lat == basedata$lat[ i ]

d2 <- subset(data, lon == basedata$lon[ j ] & lat == basedata$lat[ j ]

で、このdに該当する数字が見つからなかったときに

(d2はヒットして5だったとすると)

d + d2 とかしたい場合

そのままだと計算できないけど

length(d) + d2 だと、5が返される。

 

もっとも、dに数字がヒットした時に困るから

その前にifとかでちゃんと仕分けしておかないといけないけど。

 

強引かもしれないけど、なんとか解決。

後方流跡線をRで描いてみる・終

後方流跡線を描くとき

つまづいたのが、geom_line と geom_pathの違い

 

ずっとlineで描いていたけど、ある線が、おかしくなった。

(例のデータがどれだったのか、解らないので

 例図がありません。そのうち見つけたら、更新)

 

なんだか、ギザギザになって明らかにおかしい。

 

で、細かく見てみると

lonのx軸が小さいほうから順番に連なっていることが発覚。

後方流跡線にありがちな、くるっと回っていくのが

うまく表現できていませんでした。

 

その解決のための変更がgeom_path

線を引くときは、どちらがいいのかちゃんと考えないと。

 

知ってる人は知っているんだろうけど

一瞬、絶望的な気分になった。

 

 

後方流跡線をRで描いてみる・続続

後方流跡線のデータを取り出す方法も分かった。

地図も描いた。

 

なので、あとはこれを

くっつけるだけ。

 

いろいろと捗るといいなぁと思って

データのファイル名を引数に関数にしてみた。

ちなみに、後方流跡線はNOAAのものです。

気象データはGDASの0.5度グリッド。

 

backtraje <- function(data) {
i <- 11
lat <- na.omit(as.numeric(unlist(strsplit(data[i]," ")))/1)[10]
lon <- na.omit(as.numeric(unlist(strsplit(data[i]," ")))/1)[11]
height <- na.omit(as.numeric(unlist(strsplit(data[i]," ")))/1)[12]

data3 <- as.numeric(c(lat, lon , height))

for (i in 12:length(data)){
lat <- na.omit(as.numeric(unlist(strsplit(data[i]," ")))/1)[10]
lon <- na.omit(as.numeric(unlist(strsplit(data[i]," ")))/1)[11]
height <- na.omit(as.numeric(unlist(strsplit(data[i]," ")))/1)[12]

data3 <- rbind(data3, as.numeric(c(lat, lon, height)))
}

colnames(data3) <- c("lat", "lon", "height")
data3 <- as.data.frame(data3)

#流跡線かくときはlineじゃなくてpath
ggmap3 <-
ggmap2 +
geom_path(data = data3, aes(x=lon, y = lat), size=0.5) +
geom_point(data = data3, aes(x=lon, y = lat, size=height), shape = 1, colour="red") +
labs(size="Height(m)") +
scale_radius(range=c(0,10),limits=c(0,max(data3$height))) +
xlab("") +
theme(axis.text.x = element_blank(),
legend.position = "top")

 

}

 

で、データを読み込む(データ名:12110103.txt)

そして関数で描画

 

data <- readLines("12110103.txt")
backtraje(data)

 

f:id:Q21:20170714182034p:plain

なんでGMTじゃないかというと

こうやって、高度の情報を組み込みたかったから。

GMTでも、できるのかもしれないけど・・・・)

 

高度が高いところは円が大きく、低いところは小さい。

上空から見た感覚でイメージしてみた。

色で分けたのも作ったけど、解りにくくなったので

大きさで高度を表現。

 

これでサクサク見ていけば

なにか感じることが出てくるかもね~

後方流跡線をRで描いてみる・続

後方流跡線のデータを読み込んだら

あとはこっちのもの(なにが?)

 

地図を描いて、後方流跡線を乗っけるだけ。

 

今回、地図はNature Earthのものを使用しました。

荒めの国境データですが、今回にはこれで十分かと。

shapeファイルを読み込んで、ggplot2で描けるように

変換しました。

 

library(ggplot2)

library(maptools)

worldCountry10_shp <-  readShapeLines("ne_10m_admin_0_countries.shp")
worldCountry10_map <- fortify(worldCountry10_shp)

ggmap1 <-
  ggplot() +
  geom_path(data = worldCountry10_map, aes(x = long, y = lat, group = group),                                    colour="#999999") +
  xlim(c(80, 150)) + ylim(c(20, 60)) +
  coord_fixed(ratio=1.35/1)

 

coord_fixedは目視で設定した、個人的な趣味の数字です。

正解は、何なんだろう?

 

あと、中国は土地勘が無いので適当な都市を地図で表示。

 

pointdata <- read.csv("pointdata.csv", header=T)

> pointdata
name lat lon
1 Beijing 39.90449 116.3915
2 Shanghai 31.24787 121.4727
3 Hong_Kong 22.27838 114.1743
4 Ulan_Bator 47.92138 106.9055
5 Shenyang 41.80030 123.4338
6 Chongqing 29.57000 106.5800

 

ggmap2 <-
  ggmap1 +
  geom_point(data = pointdata,aes(x = lon, y = lat), size=4, shape=16, colour="blue") +
  geom_text(data = pointdata, aes(x= lon+1, y = lat-1, label = name), size=5)

 

f:id:Q21:20170714101958p:plain

この段階で地図を描くと、こんな感じ。

 

あとは、これに後方流跡線のデータを載せるだけ。

 それはまた、次に。

後方流跡線をRで描いてみる

オフラインで得た後方流跡線のテキストデータ。

これまではExcelマクロ&GMTで描いていたけど

様々な処理をRにまとめたくて

Rで描くことに挑戦。

 

で、まず最初に戸惑ったのがデータの読み込み。

今まではExcelで読み込んでいたから問題なかったけど

Rで行ごとに項目数が異なるデータを読み込むことが

簡単ではなかったので、その覚書。

 

つまりは、read.tabel()では、読み込めませんでした。
 
なので、readLines()で読み込みます。
ただし、これだと1行が一つの文字列として読み込まれるので
それぞれの数字を、ばらばらにしないといけない。
 
そこで使ったのがunlist(strsplit(data[i]," "))
空白を区切りとして、それぞれの行を個別の文字に認識。
 
しかし、ここでまた問題が。
なんと、各文字間の空白の数が、ばらばらだったのです。
つまり、
1□□□□□1
1□□□□ - 1
1□□□ 1 2 5
 
のように、後ろの数が大きくなると間の空白の数が変わるのです。
(なんて気の利いた設定だこと・・・・・・)
 
そうなると、空白を区切り文字にすると
各要素へのアクセスをするときに
ある行は、目的の数字がdata[48]だったとしても
別の行では、同じ緯度なりがdata[59]だったりするわけで
目的の要素にアクセスできるとは、限らなくなるわけで。
行ごとに、指定する数が変わってきてしまう!
 
なんとか空白の要素だけを取り除きたいけどis.na()とかを使っても
空白はチェックを通ってしまう。
 
なので一工夫。
要素を数字に変換して、1で割ってやると
数字はそのまま、空白はNAになるのですよ!!
で、このNAを除けば、数字だけが残るという仕組み。
 
つまり、後方流跡線のファイルをsampldata.txtとすると
 
data <- readLines("sampledata.txt")
i <- 11 #これは一例。目的の行数を入れる
lat <- na.omit(as.numeric(unlist(strsplit(data[11]," ")))/1)[9]
lon <-na.omit(as.numeric(unlist(strsplit(data[11]," ")))/1)[10]
height <- na.omit(as.numeric(unlist(strsplit(data[11]," ")))/1)[11]
 
data2 <- c(lat, lon, height)
 
これで、緯度・経度・高度のベクトルが得られるね!