2024-04-20

メソ数値予報モデル(MSM)のデータ領域可視化(その1)

Python Debian 気象学

メソ数値予報モデル(MSM)のデータ領域を可視化してみました。まずは結果から貼り付けておきます。

必要なパラメータさえ分かればMSMのデータ領域は自分で計算することはできますが、ここでは、公開されているMSMサンプルデータ(GRIB2)をpygribを用いて読み込み、緯度経度の情報を取り出して図示しています。必要なものは以下のとおりです。

MSMサンプルデータ(GRIB2) 気象庁 | GPVサンプルデータの一覧
Pythonライブラリ(pygrib, numpy, matplotlib)

Pythonライブラリについては、以下のimportが成功すれば準備OKです。

import pygrib
import numpy as np
import matplotlib.pyplot as plt

Debian/Ubuntu系であれば以下のようにaptでインストールできるでしょう。

sudo apt install -y python3-grib python3-numpy python3-matplotlib

コードは以下のとおりです。極めて単純です。気をつけないといけないのは、経度方向にnx, 緯度方向にny格子点があるとき、latlons() の返り値 lats, lons は (ny, nx) の2次元配列になっていることですかね(慣れない)。

# -*- coding: utf-8 -*-
# MSMサンプルデータをpygribで読み込み、データ領域を描画する
import pygrib
import numpy as np
import matplotlib.pyplot as plt

def extract_latlons(filename):
    data = pygrib.open(filename)
    lats,lons = data.message(1).latlons() # 読み込む物理量はなんでもよい
    nlats = np.concatenate([lats[0,:],lats[:,-1],lats[-1,::-1],lats[::-1,0]])
    nlons = np.concatenate([lons[0,:],lons[:,-1],lons[-1,::-1],lons[::-1,0]])
    return nlats,nlons

if __name__ == '__main__':
    file1 = "Z__C_RJTD_20240228000000_MSM_GPV_Rjp_Glm5km_Lm1-39_Pqc_FH00-39_grib2.bin"
    file2 = "Z__C_RJTD_20171205000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin"
    lats1,lons1 = extract_latlons(file1)
    lats2,lons2 = extract_latlons(file2)
    
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.set_aspect("equal")
    ax.set_xlim(100,160)
    ax.set_ylim(15,60)

    ax.plot(lons1,lats1,color="red",label="Lambert (817x661)")
    ax.plot(lons2,lats2,color="blue",linestyle=":",label="LatLon (481x505)")

    ax.text(0.01,-0.1,__file__,color="blue",fontsize="smaller",
            transform=ax.transAxes)
    
    plt.legend()
    plt.title("Meso-Scale Model's domain")
    plt.savefig("figure.png",dpi=150)
    plt.close

なお、今回の図示は簡素化したものであって、海岸線とかを一切描画していません。続編として今度はCartopyで凝った図を描いてみたいと思いますが今回はこのへんで。

2024-01-20

ScanSnap iX1500を使用してPC2台でのスキャン比較結果からカラー濃度の差異が判明

ScanSnap

　ScanSnap iX1500を使用している者です。もともとNEC Lavieでスキャンをしていましたが、NEC Lavieの動作が遅くなってきて使用に耐えられなくなったので、もう1台所持しているLenovo PCをメインに使用することにしました。同一のScanSnapでPCを複数台使うためにはScanSnap Cloudにユーザー登録の必要があり、かつ2台目以降、スキャンのプロファイルの設定をしないといけなかったので、ユーザー登録をして、ここ半年以上使っていました。
　ところが、ここ数ヶ月になってLenovo PCでカラースキャンした結果の色の濃さが薄くなる症状に悩まされるようになりました。あるときまでは特に違和感なくスキャンできていたと感じていたので、それが変わった原因が分からず、しばらく悩み、本来の消耗品の寿命かとも考えましたが、原因を究明する方法は、NEC LavieとLenovo PCの2台で同一の書類をスキャンして比較するしかないと思うに至りました。
　古いPCの起動に少々手こずったものの、同一の書類をスキャンした結果、カラー濃度に差が見られた一方で、ScanSnapにデフォルトで搭載されているプロファイル設定でスキャンすると差異がほとんど見られなかったことから、PC2台のプロファイルの設定に微妙な差があると分かりました。その結果、「裏写りを軽減します」にチェックが付いていたことがカラー濃度を薄くしていた要因と分かりました。
　ということでまとめると、カラー濃度が薄いと思った際には「裏写りを軽減します」にチェックが付いてしまっているかどうかを確認いただければ改善する可能性があります。お悩みの方はお試しあれ。
　今回の教訓としては、何かバージョン等が変わったときに挙動が変化していないかどうかを確認するために、ベンチマークとなる指標を用意しておくべきだと思いました。

2023-12-17

Windows11 PC の WSL2環境から Intel OneAPI をアンインストールする

ifort

以前の記事で Intel Fortran を利用するという記事を書きました。
matsubi0526.hatenablog.com
記事の内容のとおり、手許にある Let's Note (Debian) と Lenovo PC (WSL2) に Intel OneAPI をインストールしたのですが、後者の Windows11 PC において、ディスク容量が非常に逼迫している状況がここ最近続いており、今回、Intel Fortran の利用頻度を鑑みてアンインストールすることにしました。アンインストールと言っても手順は簡単で、インストールしたパッケージを sudo apt autoremove packagesすることでできました。ディスク容量をこれで約17GBも空けることができました。ディスク領域使用しすぎですね。とりあえず intel-basekit を指定するだけで自分の場合、同時に intel-hpckit もアンインストールされました。

sudo apt autoremove intel-basekit

2023-12-16

高校への数学12月号(高数オリンピック)

数学 Ruby

(n-1)人目までの得票結果に対して、n人目の得票結果を反映させる、というのがミソです。これを一気にn人の得票パターンを網羅的に計算させようものなら、n=8くらいからPCのメモリが不足してしまい強制終了になります。実際、自分は最初にそのように素朴なコードを組んで検算していました。方針を転換して一度外部CSVに吐き出し、bash側でuniq処理をかけて再度Rubyで読み込むという操作をしてようやくn=9のときの計算結果を出すことができましたが、実行には150分ほどかかってしまいました。帰納的に処理するやり方では1秒もかかりません。やり方の工夫次第では計算量を大幅に削減できるものだと気付かされる一例になりました。
なお、解析的に計算すると、(1)は ${}_{2n+3}C_{3} = \dfrac{1}{3}(n+1)(2n+1)(2n+3)$ 、(2)は ${}_{2n+3}C_{3} - 4{}_{n+2}C_{3} = \dfrac{1}{3}(n+1)(2n^{2}+4n+3)$ となり、Ruby スクリプト中では解析解との比較を実施しています。

# -*- coding: utf-8 -*-
module Define
  module_function
  def make_array(ylist)
    xlist = Array.new
    xlist[0] = ylist.count("A")
    xlist[1] = ylist.count("B")
    xlist[2] = ylist.count("C")
    xlist[3] = ylist.count("D")
    return xlist
  end
  # 決め打ちで9回計算を前提とする
  def calc_answer(vlist)
    clist = Array.new
    rlist = Array.new
    for aa in vlist
      xx = Define::make_array(aa)
      rlist.push(xx)
    end
    rlist.uniq!
    clist[0] = rlist.length
    
    # 2回目以降は帰納的に計算できる
    for ii in 2..9
      ylist = Array.new
      for bb in vlist
        xx = Define::make_array(bb)
        for rr in rlist
          yy = [rr,xx].transpose.map{|ary| ary.inject(:+)}
          ylist.push(yy)
        end
        ylist.uniq!
      end
      clist[ii-1] = ylist.length
      rlist = ylist
    end
    return clist
  end
end

if __FILE__ == $0 then
  # --- Question 1 --- #
  vlist1 = ["A","B","C","D"].repeated_combination(2).to_a
  p Define::calc_answer(vlist1)
  p (1..9).to_a.map{|ii| (ii+1)*(2*ii+1)*(2*ii+3)/3}

  # --- Question 2 --- #
  vlist2 = ["A","B","C","D"].combination(2).to_a
  p Define::calc_answer(vlist2)
  p (1..9).to_a.map{|ii| (ii+1)*(2*ii*(ii+2)+3)/3}
end

途中でやっている Array#transpose.map{|ary| ary.inject(:+)} は配列の要素どうしの和を取りたいがための操作。
qiita.com

2023-12-03

石橋を叩いたら

「石橋を叩いて結局渡らない」と揶揄されることがあったので、やらなかった後悔よりもやった後悔を大事にしたい。

2023-11-17

プログラムのコードを書くならQiitaのほうがいいのか

コードを掲載することが多々あるけれども、それならばQiitaに乗り移ったほうがいいのかしら。どうなんでしょう。

2023-09-14

findで探索した結果に対してchmodで権限付与する

bash

パーミッションが444(-r--r--r--)なファイルを検索して、write権限を付ける(-rw-r--r--)ときのコマンド。findでパーミッションが444のファイルを探索し、それらに対してchmodを実行するということをやっています。

$ find -type f -perm 444 | xargs chmod 644

2023-09-10

二番煎じな情報が巷には多いことについて

似たような記事は淘汰されていくような仕組みがないものだろうか。

2023-08-13

文字列内の複数パターン置換をHashで行う

Ruby

半年くらい前に自作したスクリプトに組み込んだ処理なのに、どこのサイトでやり方を調べたのか忘れてしまった。たぶん、「Ruby 置換 Hash」とかで調べて、なるほどふむふむと思いながら実装したんだろう、String#gsub(pattern, hash) とすることで複数置換ができる。

str = "first second third"
gsubhash = { "first"=>"Mike","second"=>"Johnny","third"=>"Tom" }
p str.gsub(Regexp.union(gsubhash.keys), gsubhash)
 #=> "Mike Johnny Tom"

gsubで置換するところ、正規表現を与えるのに Regexp#union とすればよいのだと、調べる過程で初めて知った。自分が当初書いていたスクリプトでは、Array#join を使ってなんか無理矢理感を出していた。

p str.gsub(/#{gsubhash.keys.join('|')}/, gsubhash)

weblog.2410.dev

2023-05-24

2023年6月号大問5

数学 Ruby

証明問題なので、この記事でやることはただの実験に過ぎません。
$a_{1} = 1, a_{2} = 2, a_{3} = 3, \;\; a_{n+3} = a_{n+2} + a_{n} \;\;(n=1,2,3, \cdots)$
によって数列 $a_{n}$ を定めるとき次式が成り立つことを示せ、というものです。
$\displaystyle a_{n} = \sum_{k=0}^{\lfloor \frac{n+1}{3} \rfloor} {}_{n+1-2k}C_{k}$
さて実験してみましょう。二項係数の計算部分については以下のサイトを参考にしました。
tsujimotter.hatenablog.com

# -*- coding: utf-8 -*-
module Define
  module_function
  # 2項係数を計算する
  def binom(n,k)
   k = [k, n-k].min
   if (k==0) then
      val = 1
   else
      val = Define::binom(n-1,k-1)*n/k
   end
   return val
  end
end
if __FILE__ == $0 then
  nmax = 10
  aa = Array.new(nmax+1)
  bb = Array.new(nmax+1)

  # --- 漸化式による計算 --- #
  aa[1] = 1; aa[2] = 2; aa[3] = 3
  for nn in 1..nmax-3
    aa[nn+3] = aa[nn+2] + aa[nn]
  end

  # --- 二項係数による総和 --- #
  for nn in 1..nmax
    sum = 0
    for kk in 0..((nn+1)/3)
      sum += Define::binom(nn+1-2*kk,kk)
    end
    bb[nn] = sum
  end
  
  # --- 結果の出力 --- #
  p aa
  p bb
  p aa==bb
end

実行してみた結果が以下のとおりです。ちゃんと成立していますね。

$ ruby main.rb 
[nil, 1, 2, 3, 4, 6, 9, 13, 19, 28, 41]
[nil, 1, 2, 3, 4, 6, 9, 13, 19, 28, 41]
true

Rubyの場合、配列の要素番号は0始まりなので、最初の項は nil になっています。
ガウス記号のところはわざわざメソッド化しなくても、整数どうしの除算で整数部分だけ取り出せるので特に工夫していません。

2023-05-10

2023年5月号学コン大問2

数学 Ruby

手計算でやると数え上げが面倒で、漏れがあるか心配だったので、Rubyで検算用のスクリプトを書きました。コードは別に難しくなく、答えは一瞬で出ます。使うメソッドとしては Array#permutation くらい。あとはテクニックとして 5桁の数字を構成するときにいわゆる Horner法を用いていることくらいでしょうか。nn を計算しているところ、工夫すれば再帰的に書けるような気がしないでもないです(桁数が上がっても対応できるような汎用性があると有益なんだろうけど)。

# -*- coding: utf-8 -*-
# 2023年5月号学コン大問2
# 無作為に抽出した5枚のカードで数字Nを作り
# Nが6の倍数になるときの確率を求める
#
if __FILE__ == $0 then
  aa = [1,2,2,3,3,3,4,4,4,4,5,5,5,5,5]
  ii = 0; mm = 0
  
  for cc in aa.permutation(5).to_a
    mm += 1
    if (cc.length == 5) then
      nn = ((((cc[0]*10)+cc[1])*10+cc[2])*10+cc[3])*10+cc[4]
      if (nn.modulo(6) == 0) then
        ii += 1
      end
    else
      raise "ary length must be 5: cc.length = #{cc.length}"
    end
  end
  p Rational(ii,mm) #=> (58/429)
end

2023-04-15

シェルスクリプトで実行しているスクリプトのパスを取得する

bash

やっぱりどこかの記事の二番煎じになってしまうけど、よく忘れるのでメモ。

ROOTDIR=$(cd $(dirname $0) && pwd)

常にファイルが置かれているディレクトリが変数 ROOTDIR に格納される。

2023-04-13

a以上b以下の連続する自然数の総和が500になるもの

数学

こちらの動画にある数学の問題。Rubyでプログラムを書いて強引に解くとどうなるか検証しました。プログラムも数行で簡単にかけてしまいます。
www.youtube.com

# -*- coding: utf-8 -*-
# 2023-04-12 created 
# https://www.youtube.com/watch?v=ScGsag-RAsA&t=357s
if __FILE__ == $0 then
  for bb in 2..500
    for aa in 1..bb-1
      sum = (aa..bb).to_a.sum
      if (sum == 500) then
        puts "(a,b) = (#{aa},#{bb})"
      end
    end
  end
end

実行結果はこちら。1秒もかからずに瞬殺できます。

$ time ruby sum500.rb
(a,b) = (8,32)
(a,b) = (59,66)
(a,b) = (98,102)

real    0m0.699s
user    0m0.688s
sys     0m0.011s

動画にあるように、論理的に絞り込んでいく思考力の方が大事ですが、検算の手法を知っておくことも大事かと思います。

2023-04-12

同値な条件式

数学

$a,b,c$ の少なくともひとつは1である、という条件を同値な数式で書けば、
$(a-1)(b-1)(c-1) = 0$
$a,b,c$ のすべてが1である、という条件を同値な数式で書けば、
$(a-1)^2+(b-1)^2+(c-1)^2 = 0$
以下の動画を見たときになるほどと思いました。自分だったら解と係数の関係から3次方程式を立式して解くかなと思いましたけど、同値な条件式を使えば楽にできますね。
www.youtube.com

2023-04-10

2023年4月号学コン大問1

Python

2通りの図示ができることがトラップ。

図示に使ったPython スクリプトも掲載しておきます。洗練された書き方にはなっていないだろうけれども、参考になればと思います。ちなみに動作環境は Python 3.9.2 です。

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as pat

def calc_gaisetsu_circle_hankei(a,b,c):
    menseki = calc_menseki(a,b,c)
    return a*b*c/(4*menseki)

def calc_menseki(a,b,c):
    s = (a+b+c)/2
    return np.sqrt(s*(s-a)*(s-b)*(s-c))

a = 5; b = 7; c = 6
CCR = calc_gaisetsu_circle_hankei(a,b,c)
print("----- check 1 -------")
print("外接円半径(数値計算) = {}".format(CCR))
print("外接円半径(解析解)   = {}".format(35/(4*np.sqrt(6))))

cosA = (b**2+c**2-a**2)/(2*b*c)
cosB = (c**2+a**2-b**2)/(2*c*a)
cosC = (a**2+b**2-c**2)/(2*a*b)
OF = np.sqrt(CCR**2-25/4)
vb = np.array((-5/2,-OF))
vc = np.array((5/2,-OF))
va = -(b*cosB*vb+c*cosC*vc)/(a*cosA)

k = 35/16 # 相似比
cp = k*c
ap = k*a
bp = k*b
sp = (ap+bp+cp)/2

print("----- check 2 -------")
print("A'B'(数値計算) = {}".format(cp))
print("A'B'(解析解)   = {}".format(105/8))

print(sp-ap)
print(sp-bp)
print(sp-cp)

vbp = np.array((-(sp-bp),-CCR))
vcp = np.array((sp-cp,-CCR))
vap = -(bp*vbp+cp*vcp)/ap

# --- Main Routine (drawing) --- #
fig = plt.figure(figsize=(12,6))
ax = fig.add_subplot(121)
ax.set_xlim(-7.5,7.5)
ax.set_ylim(-5,10)

circCCR = pat.Circle(xy=(0,0),radius=CCR,color="lime",ec="darkred")
p1 = pat.Polygon(xy=[va,vb,vc],fc="lightgrey",ec="darkred")
p2 = pat.Polygon(xy=[vap,vbp,vcp],fc="lightblue",ec="darkblue")
ax.add_patch(p2)
ax.add_patch(circCCR)
ax.add_patch(p1)

# --- 原点O --- #
ax.plot(0,0,marker="o",color="black")
# --- ABC --- #
ax.plot(va[0],va[1],marker="o",color="red")
ax.plot(vb[0],vb[1],marker="o",color="red")
ax.plot(vc[0],vc[1],marker="o",color="red")
# --- A'B'C' --- #
ax.plot(vap[0],vap[1],marker="o",color="blue")
ax.plot(vbp[0],vbp[1],marker="o",color="blue")
ax.plot(vcp[0],vcp[1],marker="o",color="blue")
ax.set_aspect('equal')

# ---- 点対称な図形も条件を満たす ---- #

ax2 = fig.add_subplot(122)
ax2.set_xlim(-7.5,7.5)
ax2.set_ylim(-10,5)

vap = -vap; vbp = -vbp; vcp = -vcp # 点対称
circCCR = pat.Circle(xy=(0,0),radius=CCR,color="lime",ec="darkred")
p1 = pat.Polygon(xy=[va,vb,vc],fc="lightgrey",ec="darkred")
p2 = pat.Polygon(xy=[vap,vbp,vcp],fc="lightblue",ec="darkblue")
ax2.add_patch(p2)
ax2.add_patch(circCCR)
ax2.add_patch(p1)

# --- 原点O --- #
ax2.plot(0,0,marker="o",color="black")
# --- ABC --- #
ax2.plot(va[0],va[1],marker="o",color="red")
ax2.plot(vb[0],vb[1],marker="o",color="red")
ax2.plot(vc[0],vc[1],marker="o",color="red")
# --- A'B'C' --- #
ax2.plot(vap[0],vap[1],marker="o",color="blue")
ax2.plot(vbp[0],vbp[1],marker="o",color="blue")
ax2.plot(vcp[0],vcp[1],marker="o",color="blue")
ax2.set_aspect('equal')

plt.savefig("figure.png",dpi=150)
plt.close