###pirikaについて
pirikaという言葉を初めて聞く、という人は今すぐアプリストアで「pirika」と検索してインストールしましょう。皆でご近所、故郷、旅行先のポイ捨てゴミを拾って写真を投稿しましょう。私も最近始めました。
このpirika、ゴミ拾いアプリというものですが、開発元では日本・地球上のゴミを減らすべく様々な活動をしています。詳細はHPに書いてありますが、その中で、日本各地のマイクロプラスチックゴミの情報がオープンデータとして公開されています。データはマイクロプラスチックの外観写真、IR(赤外分光)スペクトルで、2018～2020調査分まで公開されています。このIRデータをケモメトリクス手法で解析したらどうなるか興味があったので、DigitalMicrographで処理してみました。ケモメトリクスとして主成分分析(Principle Component Analysis, PCA)をするのにtemDMプラグインを使います。
###IRデータの読み込み
公開データは年度ごとに微妙にフォーマットが違います。今回は2020年のデータを対象にしました。スペクトル処理をするために、csvファイルを読み込み順番に並べる必要があります。試行錯誤しながら書いたscriptがこちら。ファイルリストをTagGroupで管理しています。
//read .csv as image
image readCSV(string path,number datatype,number header)//datatype=2 for real4
{
number fileID=OpenFileForReading(path)
object fstream=NewStreamFromFileReference(fileID,1)
number n=0
string text
while(n<header)
{
fstream.StreamReadTextLine(0,text)
n=n+1
}
number bLinesAreRows=1
number bSizeByCount=1
object imgSizeObj=Alloc("ImageData_ImageDataSize")
image img:=ImageImportTextData(path,fstream,datatype,imgSizeObj,bLinesAreRows,bSizeByCount)
return img
}
//add tag to an image
void addTg(TagGroup Tg,image img)
{
TagGroup imgTg=img.ImageGetTagGroup()
TagGroupCopyTagsFrom(imgTg,Tg)
}
//get path of target directory
string filename
string targetname="ATR.csv" //"ATR.CSV" for 2019 data (uppercase character)
if(!OpenDialog(Null, "Select "+targetname+" file", "*"+targetname, filename)) exit(0)
string dirpath=GetApplicationDirectory(0,0)
result("exploring the below directory\n")
result(dirpath+"\n")
//get file list in the target diectory
TagGroup fileList=GetFilesInDirectory(dirpath,1)
//filelist.TagGroupOpenBrowserWindow("file list",0)
//get file list of "*ATR.csv"
TagGroup targetList=NewTagList()
number n
string tempname
for(n=0;n<fileList.TagGroupCountTags();n++)
{
//check file name
fileList.TagGroupGetTagAsString("["+n+"]:Name",tempname)
if(tempname.len()>targetname.len())
{
if(tempname.right(targetname.len())==targetname)
{
targetList.TagGroupInsertTagAsString(infinity(),tempname)
}
}
}
result("number of ***"+targetname+" files: "+n+"\n")
//targetList.TagGroupOpenBrowserWindow(targetname,0)
//opening csv files if file format is the same as the first one
//opened csv image is stored into a tag group "openListAsImg".
string path //file path
number datatype=2 //datatype=2 for real4 (float)
number header=1 //number of header rows in csv
image firstImg //image of the first csv file in the targetList
TagGroup openList=NewTagList() //list of the name of the opened files
TagGroup openListAsImg=NewTagList() //image list of the opened files
number x,y,x0,y0
//first csv image
targetList.TagGroupGetIndexedTagAsString(0,path)
firstImg:=readCSV(dirpath+path,datatype,header)
firstImg.GetSize(x0,y0)
openListAsImg.TagGroupInsertTagAsArray(infinity(),firstImg)
openList.TagGroupInsertTagAsString(infinity(),path)
//opening csv files
for(n=1;n<targetList.TagGroupCountTags();n++)
{
targetList.TagGroupGetIndexedTagAsString(n,path)
//check size of the csv file
readCSV(dirpath+path,datatype,header).GetSize(x,y)
if((x==x0)&&(y==y0))
{
openListAsImg.TagGroupInsertTagAsArray(infinity(),readCSV(dirpath+path,datatype,header))
openList.TagGroupInsertTagAsString(infinity(),path)
}
}
//stack all opened csv data
number z0=openListAsImg.TagGroupCountTags()
image stackImg:=RealImage(targetname,4,x0,y0,z0)
for(n=0;n<z0;n++)
{
openListAsImg.TagGroupGetIndexedTagAsArray(n,stackImg[0,0,n,x0,y0,n+1])
}
//set openList tag to "file list" tag
TagGroup openFileList=NewTagGroup()
openFileList.TagGroupSetTagAsTagGroup("file list",openList)
//image stackImg should have (x,y,z)=(col,row,number of opened data)
result("number of opened data files: "+z0+"\n")
//wavelength as an image
image wavelength=stackImg.slice2(0,0,0,1,y0,1,2,z0,1)
wavelength.SetName("wavelenghth")
addTg(openFileList,wavelength)
wavelength.showimage()
result("wavelength range: "+wavelength.min()+" - "+wavelength.max()+"\n")
//spetrum as an image
image spc=stackImg.slice2(1,0,0,1,y0,1,2,z0,1)
spc.SetName("spectrum")
addTg(openFileList,spc)
spc.showimage()
//spectrum as an spestrum imaging type data (3D)
image spccube=stackimg.slice3(1,0,0,2,z0,1,0,1,1,1,y0,1)
spccube.SetName("spectrum cube")
addTg(openFileList,spccube)
spccube.showimage()
###データ前処理
出力した2020年のデータを見ると絶対値が1を超えるものがたまにあります。出力データの算出方法が分かりませんが、透過率とすると1を超えるのはおかしいように思います。データとしておかしくなくとも、多くのデータが1未満で、主成分解析は外れ値に弱そうなので、絶対値が1を超えるデータはいったん無視することにします。負の値が出るのもよく分からないですが、負のデータはいったんそのままにしました。
//delete data which has a value larger than one
image in:=getFrontImage() //"spectrum cube" should be the front.
number x,y,z
in.Get3DSize(x,y,z)
string name
in.GetName(name)
TagGroup inTg=in.ImageGetTagGroup()
TagGroup originalFileList
inTg.TagGroupGetTagAsTagGroup("file list",originalFileList)
TagGroup cpTg=originalFileList.TagGroupClone() //copy
//store image into tag group
TagGroup imgTg=NewTagList()
image temp
number n
for(n=x-1;n>=0;n--)
{
temp=in.slice1(n,0,0,2,z,1)
if(temp.abs().max()>1)
{
cpTg.TagGroupDeleteTagWithIndex(n)
}
else
{
imgTg.TagGroupInsertTagAsArray(0,temp)
}
}
number nn=cpTg.TagGroupCountTags()
image delImg:=RealImage(name+" del",4,nn,y,z)
for(n=0;n<nn;n++)
{
imgTg.TagGroupGetIndexedTagAsArray(n,temp)
delImg[n,0,0,n+1,y,z]=temp.slice3(0,0,0,1,1,1,2,1,1,0,z,1)
}
TagGroup resultTg=delImg.ImageGetTagGroup()
resultTg.TagGroupSetTagAsTagGroup("file list",cpTg)
delImg.showimage()
###PCA
ここで生成したデータに主成分分析をかけます。データさえ作ってしまえばtemDM使ってクリックぽちぽちするだけで解析完了です。scree plot (variance) 見るとざっくり第3成分までが主要なスペクトルですが、Loadの形を確認するとかなり高次までスペクトル情報が残っていました。ちなみに、varianceはlog表示されるようなのでこのデータでは負の値になってます。
ここで出てくる主成分の解釈はどうもよくわかりませんでした。第1主成分(PCA1)はあまりきれいなスペクトルの形をしておらずやけに斜めなので、主にバックグラウンドの強度を表す成分じゃないかと考えています。PCA2, PCA3は正負を無視すればピークの位置は似ていて、どちらもPEぽいスペクトルに見えます。因子回転すると一見PEスペクトルとバックグラウンド成分に分かれましたが、2種類のPEスペクトルと解釈すべきなのかよくわかりません。PC4はPPのスペクトルを反映しているように見えます。PC5以降はどう解釈したらよいかさっぱりわかりません。pirikaでまとめた資料を見ると、PE、PP以外にPS,Nylon等もそれなりの割合で含まれているようですが、PCAでこれらのスペクトルを抽出するのは難しそうです。データ処理まではできましたが、解釈できないのでここまでにします。ゴミを分類してやる！と意気込んでましたが、私の力では役に立ちませんでした。残念。
###PCAその2
おまけです。ここではtemDMを使ってPCAの計算をしましたが、DigitalMicrographの新しいverではGatan独自のMSA (multi-variate statistical analysis)機能が増えておりPCAも可能でした。GMS3.4以降の新しい機能のようなので、GMS3.xの初期だとたぶんできないです。temDMだと入力データを3次元にしないと動かないですが、Gatanの関数だとスペクトルを並べただけの2次元データでも入力を受け付けました。
image spc:=GetFrontImage() //2D image can be loaded
image score,load,scree
MSA_PCA(spc,score,load,scree,0)
scree.setname("scree")
score.setname("score")
load.setname("load")
scree.showimage()
score.showimage()
load.showimage()
出力結果を見ると、scoreとloadの名前が入れ替わっているので、MSA_PCA()関数に入れる順番は入れ替えたほうが良いかもしれません。Gatan公式のドキュメントではこの順番です。
void MSA_PCA( BasicImage src, BasicImage scoreStack, BasicImage loadStack, BasicImage eigenValuesLog10, Boolean
useDoublePrecision )
PCAだけでなく、varimax回転やICA (independent component analysis / blind source separation) も計算できるようなので結構使えそうですが、正直temDMあればいいかな．．．