GASとPython3で日照時間の履歴を記録する

# 1.はじめに
## 1.1 きっかけ

屋根に太陽光パネルが乗っていると日射量や日照時間が気になってくる（賃貸の我が家には乗っていないが）。おそらく太陽光パネルの発電量に一番重要な情報は日射量の情報だと思うが、実は全天日射量の計測は全国で49カ所でしか行われていない。
[全天日射量の観測地点一覧　（2017年3月時点）](http://k-ichikawa.blog.enjoy.jp/etc/HP/htm/jmaP1.html)

日射量のデータベースはNEDOのものが有名で、アメダスで観測された日照時間などから日射量への換算を行って提供されている。アメダスは全国に約800拠点あるので各地での太陽光発電のシミュレーションにはこのデータが使われることが多いらしい。
[NEDOの日射量データベース](https://www.nedo.go.jp/library/nissharyou.html)


しかし、アメダスの拠点だって自分の住んでいるところからは結構遠い場合も多い。
自分の家の真上の日射情報、Web経由で取得できないだろうか？と思ったのがきっかけ。
結論から言えば、日射量の情報取得はできていない。
ただ、任意の地点での日照時間の記録は取得できるようになった。

これまでExcel VBAでマクロを作った経験があるくらいでGASもPythonもほとんど触ったことがなかったので勉強も兼ねて挑戦した記録。
日照時間から日射量に換算するためのモデルは別の機会に挑戦してみたい。


## 1.2 日照時間のデータについて

日射量と同様、日照時間も[気象庁の過去データ](http://www.data.jma.go.jp/obd/stats/etrn/index.php)からは限られた地点でしか情報を取得できない。しかし、気象庁のHPには日照時間に関して[以下のような記述](http://www.data.jma.go.jp/obd/stats/data/mdrr/man/sampling.html#amedassun2021)がある。
>（２）2021年の「推計気象分布（日照時間）」による推計値の導入
>アメダスの日照時間について、2021年3月2日から、気象衛星観測のデータを用いた「推計気象分布（日照時間）」から得る推計値としました。 日照計による観測値と推計値は性質が異なるため、それまでの日照計による観測値との統計を切断し、平年値については推計値相当の値に補正しています。

つまり、2021年3月以降は[推計気象分布（1時間ごとの日照時間）](https://www.data.jma.go.jp/bunpu/?000&suns1h)で公開されている画像をうまく使えば、気象庁のアメダス地点で公開されているものと同じ日照時間の情報が1kmメッシュ単位で得られるということになる。


#2. 任意地点での日照時間データの取得
##2.1 環境
Windows 10
Anaconda + Python3.8

## 2.2 GoogleAppsScriptで日照時間の画像取得
まずは[気象庁の推計気象分布（1時間ごとの日照時間）](https://www.data.jma.go.jp/bunpu/?000&suns1h)から日照時間の画像情報をとってくる。
各都道府県の推計日射時間の図は以下の3つで構成されていることが分かる。

背景色・地域境界図・積算日照時間（1時間）　　※地域で東京を指定した場合
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1637818/8178bec5-ef73-f2a4-f793-8c6b71646a82.jpeg" width=15%>　<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1637818/97db7319-e51e-d391-eae4-e15a8d920e81.png" width=15%>　<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1637818/23d4498d-11d6-8fd1-5841-2ad7db56f4f5.png" width=15%> 
地域境界図と積算日照時間図はそれぞれ520×520ピクセルのpngファイルになっている。この中で必要になるのは積算日照時間の図なのでこれを取得してGoogle Driveに保存する。積算日照時間図は現在を起点に48時間前までの情報が提供されているため、1日1回実行して24時間分取り出してやればいい。

~~~Javascript
function SunshineDurationPngSave(){　　　//気象庁の推計気象分布から前日24時間分の日照時間の画像を取得してGoogleドライブ保存する

  var date = new Date(); //現在日時のDateオブジェクトを作る
  var day = date.getDate();　//現在の「日」を取得
  date.setDate(day-1);　//前日日付にしたいので-1する
  var yesterday = Utilities.formatDate(date, 'JST', 'yyyyMMdd');　  //日付の表示形式を整形

for(var i=0; i<=23; i++){ 　
    var url = 'http://www.data.jma.go.jp/bunpu/img/suns1h/****/suns1h_****_' + yesterday + Utilities.formatString("%02d", i) + '00.png';
    // "****" には表示したい地区の画像IDを入れる。例えば東京なら「3016」、高知なら「4034」など
    var options = {
      'method' : 'GET',
      'muteHttpExceptions' : true,
       }
    var response = UrlFetchApp.fetch(url, options);
    var filename = yesterday + Utilities.formatString("%02d", i) + '.png';
    var fileBlob = response.getBlob().setName(filename);
    var folder = DriveApp.getFolderById('GoogleドライブのフォルダIDを入力'); 
    folder.createFile(fileBlob);
  } 
 }
~~~


##2.2 Python3で日照時間の履歴算出


取り出したい地点の座標を指定する必要があるため、地域境界図から指定ポイントのピクセル値を読み取る（今回はペイントで愚直に読みに行った）。例えば皇居なら地域境界図の<font color="red">赤十字</font>で示したあたりになる。この地域境界図（520×520ピクセル）中、<font color="red">赤十字</font>に相当するのは「402,233」なのでこれを覚えておく。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1637818/62aaad84-1901-7b5a-ae45-641e00960d8b.png" width=20%>
一定期間の積算日照時間の図がGoogle Driveに溜まったら、ローカルにダウンロードしてきて1つのフォルダに保存し、以下を実行。これで任意の地点の日照時間の履歴が取れる。


~~~Python3
import os
import cv2
import glob

path = r'C:\Pictures\*****' #積算日照時間の画像データの入ったフォルダを指定
os.chdir(path)

f = open('SunshineDuration.txt', 'a')　#保存先のファイルを開く

DailySunshineDuration=0

for fn in glob.glob('*.png'):
    file_name = fn
    img = cv2.imread(file_name)
    pixelValue = img[***,***] #取得した画像のうち日照時間を取得したい地点のピクセル値を入力（皇居なら402,233）

    if pixelValue[0]==120:
        SunshineDuration=0
    elif pixelValue[0]==154:
        SunshineDuration=0.1
    elif pixelValue[0]==188:
        SunshineDuration=0.3
    elif pixelValue[0]==244:
        SunshineDuration=0.5
    elif pixelValue[0]==136:
        SunshineDuration=0.7
    elif pixelValue[0]==0:
        SunshineDuration=0.9
    else:
        SunshineDuration=1.0

#    BGRの凡例
#    [120 120 120] 0
#    [154 154 154] 0.1
#    [188 188 188] 0.3
#    [224 224 224] 0.5
#    [136 230 240] 0.7   
#    [  0 184 248] 0.9
#    [ 56  80 244] 1.0

    os.remove(file_name)
    f.write(fn[:4]+"-"+fn[4:6]+"-"+fn[6:8]+" "+fn[8:10]+":00"+" "+str(SunshineDuration)+'\n')

f.close()
~~~


# 感想
全国、任意の地点での日照時間の履歴を取得できるようになった。実際、太陽光発電に重要な情報は日照時間の履歴ではなく「全天日射量の予測値」なのだと思うので、取得できてどれくらい意味のあるデータなのかはちょっとよくわからない。しかし、日照時間（＋気温、湿度 etc.）から日射量の推定ができるモデルはあるようなので、任意の地点の気象情報を[OpenWeatherMap](https://openweathermap.org/)なんかで長期的に取得しておけば任意の地点の全天日射量も推定できるようになるのではないだろうか…。

# 参照
[Google Apps Scriptで簡単スクレイピング](https://haloechoes.com/system-development/easy-scraping-of-google-apps-script/)