# 2018年8月期 AITCセミナー「気象データを"R"で可視化してみよう!”用教材 # イベント詳細 http://aitc.jp/events/20180803-Seminar/info.html # ※注意事項:【ユーザー名】をログインしているユーザー名に修正 # =============================================================== # サンプルデータ「iris」を触ってみる # ■3.2.R言語の基本的な文法 # 変数の参照 iris; # サンプルデータ「iris」(全150件)を表示 data(); # 全サンプルデータの一覧表示 View(iris); # 参照窓で「iris」を参照 head(iris); # 「iris」を先頭6件だけ表示 head(iris, n=10); # 「iris」を先頭10件だけ表示 iris[iris$Species=="setosa",]; # 条件に一致する行だけ取り出す iris[5.0<=iris$Sepal.Length & iris$Sepal.Length<=5.1,]; # 複数条件 iris[c(1,2,3),]; # 先頭3行だけ取り出す iris[1:3,]; # 先頭3行だけ取り出す iris[,c(3,4,5)]; # 特定の列だけ取り出す iris[,3:5]; # 特定の列だけ取り出す iris[,c("Species")]; # 列名を指定して取り出す # 変数の操作 aaa <- iris; # サンプルデータ「iris」を変数「aaa」にコピー aaa <- iris[1:10,]; # 1〜10行だけを取り出してコピー aaa <- iris[iris$Species=="setosa",]; # setosa だけを抽出 aaa; # 確認 aaa$label <- 0; # label列を追加し、初期値を0にする head(aaa); # 確認 aaa[aaa$Sepal.Length<5,]$label <- 1; # がく片の長さが5より小さいものだけ1 head(aaa); # 確認 # ■3.3.R言語の基本的な文法 # 関数の実行 summary(iris); # 変数「iris」のサマリー表示 hist(iris$Petal.Length); # 変数「iris$Petal.Length」のヒストグラム sd(iris$Petal.Length); # 変数「iris$Petal.Length」の標準偏差(sd) ? sd # 関数sdのヘルプ table(iris$Species); # 種類をカウント # 描画してみる # plotの説明:http://cse.naro.affrc.go.jp/takezawa/r-tips/r/53.html plot(iris); # 変数「iris」を描画 plot(x=iris$Petal.Length, y=iris$Petal.Width); # 花びらだけを描画 pie(table(iris$Species)); # 種類ごとの件数での円グラフ # ■3.4.描画の情報量を増やす−1 # 色を加える plot(iris, col=iris$Species); # 種類ごとに色を塗る plot(x=iris$Petal.Length, y=iris$Petal.Width, col=iris$Species); # 花びらだけ # マーカー、もしくは文字を加える plot(x=iris$Petal.Length, y=iris$Petal.Width, col=iris$Species, pch=as.numeric(iris$Species)); # 数字を指定するとマーカー plot(x=iris$Petal.Length, y=iris$Petal.Width, col=iris$Species, pch=as.character(iris$Species)); # 文字列を指定すると、先頭の1文字 plot(x=iris$Petal.Length, y=iris$Petal.Width, col=iris$Species, pch=c("s","v","x")[as.numeric(iris$Species)]); # 改良版 # ■3.4.描画の情報量を増やす−2 # 文字の大きさを変える plot(x=iris$Petal.Length, y=iris$Petal.Width, col=iris$Species, pch=c("s","v","x")[as.numeric(iris$Species)], cex=as.numeric(iris$Species)); # 文字サイズの倍率を指定(1,2,3) # もっと複雑な描画の例 # https://momonoki2017.blogspot.com/2018/05/r008-r1.html # =============================================================== # 気象データを触ってみる # 気象庁: https://www.data.jma.go.jp/gmd/risk/obsdl/index.php # 地点:「東京」→「東京」 # 項目:「月別値」、「月平均気温」「月最高気温」「月最低気温」 # 期間:1872年1月〜2018年7月 # ファイル名変更:data.csv -> data_tokyo_month.csv # ■4.4.気象データをRで読んでみる # CSVファイルを読み込んで、扱い易くする temp <- read.table("/Users/【ユーザー名】/Downloads/data_tokyo_month.csv", skip=6, header=F, sep=",", fileEncoding="SJIS"); # temp <- read.table("http://cloud.aitc.jp/20180803_R/data_tokyo_month.csv", skip=6, header=F, sep=",", fileEncoding="SJIS"); head(temp); temp <- temp[,c(1,2,5,8)]; # 使う列だけにする head(temp); names(temp) <- c("DATE","AVGTEMP", "MAXTEMP", "MINTEMP"); # 列名を付ける head(temp); temp$DATE <- paste(temp$DATE,"1", sep="/"); # yyyy/mm -> yyyy/mm/1 head(temp); temp$DATE <- as.POSIXlt(temp$DATE, format="%Y/%m/%d"); # DATEを日付型にする temp$YYYY <- as.numeric(format(temp$DATE, "%Y")); # 年だけにする temp$MM <- as.numeric(format(temp$DATE, "%m")); # 月だけにする temp <- temp[!is.na(temp$AVGTEMP),]; # 空行を消す(無い場合もある) head(temp); # 概要を確認 summary(temp); # サマリー hist(temp$AVGTEMP); # ヒストグラム # 今までで一番熱かったのはいつか? # → 最高気温が一番高かった max(temp$MAXTEMP); # 最高気温 temp[temp$MAXTEMP==max(temp$MAXTEMP),]; # その行を表示 # ★課題:一番寒かった(min())のはいつか? # → 最低気温が一番低かった # ★★★課題:自分の出身地で、一番熱かったのはいつか? # 平均気温、最高気温、最低気温を折れ線グラフで描画 # X軸は、全データが同じなので、指定しなくてOK # Y軸は、データによって違うので、固定する # 2つ目以降を重ねて描画:par(new = TRUE) work <- temp; # work にコピー head(work); # 中身を確認 plot(x=work$DATE, y=work$AVGTEMP, type="l", col="BLACK", ylim=c(0,40)); par(new = TRUE);plot(x=work$DATE, y=work$MAXTEMP, type="l", col="RED", ylim=c(0,40)); par(new = TRUE);plot(x=work$DATE, y=work$MINTEMP, type="l", col="BLUE", ylim=c(0,40)); # 2000年以降だけを描画 work <- temp[temp$YYYY>=2000,]; # 条件を指定して抽出 head(work); # 中身を確認 plot(x=work$DATE, y=work$AVGTEMP, type="l", col="BLACK", ylim=c(0,40)); par(new = TRUE);plot(x=work$DATE, y=work$MAXTEMP, type="l", col="RED", ylim=c(0,40)); par(new = TRUE);plot(x=work$DATE, y=work$MINTEMP, type="l", col="BLUE", ylim=c(0,40)); # ★課題:Y軸を固定しないとどうなるか? # ylim=c(0,40) を削除して実行 # 7月だけを描画 # 後の描画の都合で、ここだけX軸は年 work <- temp[temp$MM==7,]; # 条件を指定して抽出 head(work); # 中身を確認 plot(x=work$YYYY, y=work$AVGTEMP, type="l", col="BLACK", ylim=c(0,40)); par(new = TRUE);plot(x=work$YYYY, y=work$MAXTEMP, type="l", col="RED", ylim=c(0,40)); par(new = TRUE);plot(x=work$YYYY, y=work$MINTEMP, type="l", col="BLUE", ylim=c(0,40)); # 年と平均気温の相関係数を出す cor(work[,c("YYYY","AVGTEMP","MAXTEMP","MINTEMP")]); # 相関係数 # 年と平均気温の回帰分析 result.lm1 <- lm(formula=AVGTEMP~YYYY, data=work); abline(result.lm1 , lwd=1 , col="BLACK"); summary(result.lm1); summary(result.lm1)$coefficients; # Estimate Std. Error t value Pr(>|t|) # (Intercept) -8.23167497 5.453166897 -1.509522 1.333874e-01 # YYYY 0.01711498 0.002800886 6.110561 9.048062e-09 # 平均気温 = -8.23167497 + 0.01711498 * 年 となる # つまり、この月の平均気温は毎年0.017度づつ上がっているということ result.lm2 <- lm(formula=MAXTEMP~YYYY, data=work); abline(result.lm2 , lwd=1 , col="RED"); summary(result.lm2)$coefficients; result.lm3 <- lm(formula=MINTEMP~YYYY, data=work); abline(result.lm3 , lwd=1 , col="BLUE"); summary(result.lm3)$coefficients; # ★課題:他の月でも試してみる # 8月 # 年と1年間の平均気温の相関 # 年(YYYY)をキーに平均値を出す # aggregateの説明: http://bioinfo-dojo.net/2017/08/15/r-aggregate/ work <- temp; m <- aggregate(list(work$AVGTEMP, work$MAXTEMP, work$MINTEMP), by=list(work$YYYY), mean, na.rm=TRUE); # 平均値(mean) head(m); names(m) <- c("YYYY","AVGTEMP","MAXTEMP","MINTEMP"); # 列名をキレイにする work <- m; head(work); work <- work[work$YYYY!=1875,]; # 一番古い1875年は6月からなので、捨てる work <- work[work$YYYY!=2018,]; # 一番新しい2018年は途中までなので、捨てる head(work); # 描画 plot(x=work$YYYY, y=work$AVGTEMP, type="l", col="BLACK", ylim=c(0,40)); par(new = TRUE);plot(x=work$YYYY, y=work$MAXTEMP, type="l", col="RED", ylim=c(0,40)); par(new = TRUE);plot(x=work$YYYY, y=work$MINTEMP, type="l", col="BLUE", ylim=c(0,40)); # 年平均と気温の相関係数と回帰分析 cor(work[,c("YYYY","AVGTEMP","MAXTEMP","MINTEMP")]); # 相関係数 result.lm1 <- lm(formula=AVGTEMP~YYYY, data=work); abline(result.lm1 , lwd=1 , col="BLACK"); summary(result.lm1); summary(result.lm1)$coefficients; result.lm2 <- lm(formula=MAXTEMP~YYYY, data=work); abline(result.lm2 , lwd=1 , col="RED"); summary(result.lm2)$coefficients; result.lm3 <- lm(formula=MINTEMP~YYYY, data=work); abline(result.lm3 , lwd=1 , col="BLUE"); summary(result.lm3)$coefficients; # ★課題:1年間の平均値(mean)ではなく、最小値(min)や最大値(max)だと、どうなる? # 特に、MINTEMPの最小値は、どうなっている? # =============================================================== # 電力データ x 気温データ # 東京電力パワーグリッド: http://www.tepco.co.jp/forecast/html/download-j.html # ダウンロードする期間:2017年 # ファイル名:juyo-2017.csv # ■6.2.電力データをRで読んでみる # CSVファイルの読み込み juyo <- read.table("/Users/【ユーザー名】/Downloads/juyo-2017.csv", skip=2, header=T, sep=",", fileEncoding="SJIS"); # juyo <- read.table("http://www.tepco.co.jp/forecast/html/images/juyo-2017.csv", skip=2, header=T, sep=",", fileEncoding="SJIS"); head(juyo); # 列名に日本語が入っていると色々と面倒なので、列名を"KW"に変更 names(juyo) <- c("DATE","TIME","KW"); head(juyo); # 概要を確認 summary(juyo); # サマリー hist(juyo$KW); # ヒストグラム juyo[juyo$KW==max(juyo$KW),]; # 最大値のあった行 # 気象データは年月日と時刻が分かれているので、まとめる juyo$DATETIME <- as.POSIXlt(paste(juyo$DATE, juyo$TIME, sep=" "), format="%Y/%m/%d %H:%M"); # DATE+TIMEを日時型にする juyo$DATETIME <- juyo$DATETIME + 0; # 型を合わせる juyo$DAY <- substring(weekdays(juyo$DATETIME),1,1); # 曜日列を追加 head(juyo); table(juyo$DAY); # 1年間の電力量の推移を描画 plot(juyo$KW, type="l"); # 全部を折れ線グラフで描画 plot(juyo[juyo$TIME=="12:00",]$KW, type="l"); # 12:00のみ描画 # ★課題:他の時間帯を描画してみる work <- juyo[juyo$TIME=="12:00",]; # 12:00だけにする work <- work[work$DAY!="土" & work$DAY!="日",]; # 土日を除外 plot(work$KW, type="l"); # 12:00の平日だけを描画 # 異常?な値になっているのは何月か? work$MONTH <- as.numeric(format(work$DATETIME,"%m")); # 月の列を追加 head(work); plot(work$KW, type="b",pch=work$MONTH); # 月が数値なので、文字列化する work$MONTH <- as.character(work$MONTH); plot(work$KW, type="b",pch=work$MONTH); # 10月以降が全部"1"なので、16進数表記にする work[work$MONTH=="10",]$MONTH <- "a"; # 10月を"a"に work[work$MONTH=="11",]$MONTH <- "b"; # 11月を"b"に work[work$MONTH=="12",]$MONTH <- "c"; # 12月を"c"に plot(work$KW, type="b",pch=work$MONTH); # もっとスマートな書き方 plot(x=work$DATETIME, y=work$KW, type="l"); # ★課題:夜間は土日の影響を受けるのか? # 気象庁: https://www.data.jma.go.jp/gmd/risk/obsdl/index.php # 地点:「東京」→「東京」 # 項目:「時別値」、「気温」 # 期間:2017年1月1日〜12月31日 # ファイル名変更:data.csv -> data_tokyo_hour_2017.csv # ■7.2.気象データをRで読んでみる # CSVファイルを読み込む temp <- read.table("/Users/【ユーザー名】/Downloads/data_tokyo_hour_2017.csv", skip=5, header=F, sep=",", fileEncoding="SJIS"); # temp <- read.table("http://cloud.aitc.jp/20180803_R/data_tokyo_hour_2017.csv", skip=5, header=F, sep=",", fileEncoding="SJIS"); head(temp); temp <- temp[,c(1,2)]; # 使う列だけにする names(temp) <- c("DATETIME","TEMP"); # 列名を付ける head(temp); # 概要を確認 summary(temp); # サマリー hist(temp$TEMP); # ヒストグラム temp[temp$TEMP==max(temp$TEMP),]; # 最大値のあった行 # DATETIMEを日時型に変換する head(temp); temp$DATETIME <- as.POSIXlt(temp$DATETIME, format="%Y/%m/%d %H:%M:%S"); # DATETIMEを日時型にする temp$DATETIME <- temp$DATETIME - 3600; # 1時間(3600秒)ずらして、電力データ側の表現に合わせる head(temp); # 1年間の気温の推移を描画する plot(temp$TEMP, type="l"); plot(temp[format(temp$DATETIME, "%H:%M")=="12:00",]$TEMP, type="l"); # 12:00(元データでは13:00)だけ # ★課題:他の時間帯を表示してみる # 13:00, 01:00, 06:00 # 1日で一番温度の高い時刻を知りたい work <- temp; # コピー work$DATE <- format(temp$DATETIME, "%Y/%m/%d"); # 日付だけの列を作る head(work); m <- aggregate(work$TEMP, by=list(work$DATE), max); # 日付で最大値を見つける names(m) <- c("DATE", "MAXTEMP"); # 列名を付ける head(m); work <- merge(work, m, by=c("DATE"), all=T); # 元データにマージ m <- work[work$TEMP==work$MAXTEMP,]; # 気温と最高気温が同じ列だけを取得 head(m); m$HOUR <- as.numeric(format(m$DATETIME, "%H")); # 時刻を描画し易いよう数値化 table(m$HOUR); # 時刻のカウント plot(x=m$DATETIME, y=m$HOUR, col="RED"); # プロットしてみる hist(m$HOUR); # ヒストグラム hist(m$HOUR, breaks=seq(0, 23, 1)); # 1時間ごとのヒストグラム # ===電力データと気温データを結合=== # ■8.1.電力 x 気象の結合−3 head(juyo); head(temp); juyotemp <- merge(juyo, temp, by=c("DATETIME"), all=T); # 横方向に結合 head(juyotemp); # 結合に失敗した行を確認 juyotemp[is.na(juyotemp$KW),]; juyotemp[is.na(juyotemp$TEMP),]; # 結合に失敗した行を削除 juyotemp <- juyotemp[!is.na(juyotemp$KW),]; # 空行を消す(無い場合もある) juyotemp <- juyotemp[!is.na(juyotemp$TEMP),]; # 空行を消す(無い場合もある) head(juyotemp); # 気温と電力データを折れ線グラフに重ねる # Y軸は自動にまかせる plot(juyotemp$TEMP, type="l", col="RED"); # 気温が赤 par(new = TRUE); plot(juyotemp$KW, type="l"); # 電力が黒 # 特定の時刻だけ描画する work <- juyotemp[juyotemp$TIME=="12:00",]; head(work); plot(work$TEMP, type="l", col="RED"); par(new = TRUE); plot(work$KW, type="l"); # ★課題:他の時間帯も描画してみる # 1:00, 7:00, 17:00 # 特定の時刻から土と日を除外する work <- juyotemp[juyotemp$TIME=="12:00",]; work <- work[work$DAY!="土",]; work <- work[work$DAY!="日",]; head(work); plot(work$TEMP, type="l", col="RED"); par(new = TRUE); plot(work$KW, type="l"); # "12:00"だけ、X軸に気温、Y軸に電力量にして描画する work <- juyotemp[juyotemp$TIME=="12:00",]; head(work); plot(x=work$TEMP, y=work$KW); # 平日を黒、土曜日を青、日曜日を赤にして juyotemp$COL <- "BLACK"; juyotemp[juyotemp$DAY=="土",]$COL <- "BLUE"; juyotemp[juyotemp$DAY=="日",]$COL <- "RED"; head(juyotemp); # 再び、"12:00"だけを描画する work <- juyotemp[juyotemp$TIME=="12:00",]; plot(x=work$TEMP, y=work$KW, col=work$COL); # 黒:平日、青:土、赤:日 # ★課題:他の時間帯も描画してみる # 1:00, 7:00, 17:00 # "12:00"の平均気温が20度以上の部分だけを抽出する work1 <- juyotemp[juyotemp$TIME=="12:00" & juyotemp$TEMP>=20,]; head(work1); plot(x=work1$TEMP, y=work1$KW, col=work1$COL); # 描画 # 平日のみに絞り込む work1 <- work1[work1$DAY!="土" & work1$DAY!="日",]; # 土日を除外 plot(x=work1$TEMP, y=work1$KW, col=work1$COL); # 描画 # 妙に値の小さいデータを確認 # (どこかから休日テーブルを持ってきた方が確実) work1[work1$KW<=3000,]; # DATETIME DATE TIME KW day TEMP COL # 2917 2017-05-02 12:00:00 2017/5/2 12:00 2919 火 22.0 BLACK <- GW # 2941 2017-05-03 12:00:00 2017/5/3 12:00 2574 水 20.9 BLACK <- GW # 2965 2017-05-04 12:00:00 2017/5/4 12:00 2538 木 22.7 BLACK <- GW # 2989 2017-05-05 12:00:00 2017/5/5 12:00 2563 金 23.6 BLACK <- GW # 7357 2017-11-03 12:00:00 2017/11/3 12:00 2917 金 21.2 BLACK <- 分化の日 # 休日のデータを削除する work1 <- work1[work1$DATE!="2017/5/2",]; work1 <- work1[work1$DATE!="2017/5/3",]; work1 <- work1[work1$DATE!="2017/5/4",]; work1 <- work1[work1$DATE!="2017/5/5",]; work1 <- work1[work1$DATE!="2017/11/3",]; plot(x=work1$TEMP, y=work1$KW, col=work1$COL); # 再描画 # 相関係数(cor)を求める cor(work1[,c("KW","TEMP")]); # KW TEMP # KW 1.0000000 0.8791156 # TEMP 0.8791156 1.0000000 # 回帰分析を行う result.lm1 <- lm(formula=KW~TEMP, data=work1); abline(result.lm1 , lwd=1 , col="BLACK"); summary(result.lm1)$coefficients; # Coefficients: # Estimate Std. Error t value Pr(>|t|) # (Intercept) 375.508 175.424 2.141 0.0344 * # TEMP 129.380 6.513 19.866 <2e-16 *** # つまり、KW = 375.508 + 129.380 * TEMP; # ★課題:明日の12:00の気温が30度の場合の、電力予測 # ★課題:明日の12:00の気温が35度の場合の、電力予測 # ★課題:明日の12:00の気温が40度の場合の、電力予測 # 元のグラフに重ねて描画 plot(x=work$TEMP, y=work$KW, col=work$COL); abline(result.lm1 , lwd=1 , col="BLACK"); # ★課題:閾値を20度→25度にすると、どうなる? # 相関係数(0.8791156)は変わるか? # 回帰の線の傾き(129.380)は変わるか? # ★課題:寒い時期(15度以下)に対してもやってみる # ★課題:以下の温度の時の、電力量を予測する # 30, 35, 40 # Rで計算する方法 round(375.508 + 129.380 * c(30, 35, 40), 1); # 計算結果を四捨五入