2019年2月6日

使用imagej Macro客製化imagej的LUT(lookup table)

這兩天在我的imagej教學影片《05 imagej的LUT》上,有人留言詢問了問題:

你好,我想请问一下如何只是显示比如说25到150之间的lut,而小于25和大于150点不显示,而且在添加calibration bar的时候也是显示25-150呢?谢谢啦!
我现在的问题是我把8bit的图片的lut设置成了jet,也就是最低的灰度值为蓝色,最高的255为红色。但是呢我的图片最高的灰度值为150,最低的为25,我想在这个范围内设置为jet的lut可以吗?也就是25的灰度值为蓝色,150的灰度值为红色。且中间的颜色变化是像jet中连续的。我的问题有点儿长,不知道是不是说明白了。谢谢您!


lut是用來視覺化圖片的工具,作用就是把8bit灰階圖片著色,lut提供的就是著色的規則表。例如。例如下圖左下角的LUT規定了「72這個灰階值,要用R0 G159 B255的顏色來著色」


不同的顏色規則就會形成不同的規則表,在imagej裡頭預設就有數十種LUT可以使用。利用Image/Color/Display LUTs 可以顯示這些LUT


其實我並不太在youtube回答留言者的問題,因為光是我Gmail裡的問題我就回答不完了,再者是youtube的留言提問,並沒有辦法附圖去說明。不過因為留言者的問題剛好跟我過去想解決的問題有關,他寫出的問題有我需要的關鍵字,所以我就投入去找資料了。

留言回應者提到的LUT名稱叫做「jet」,我翻了好久倒是沒看到。後來去查才知道是Matlab這個程式裡頭的LUT,它的顏色配置就是本文第一張圖左上角的樣子,就像是彩虹的顏色。


要解決留言回應者的問題前,得先找到方法去做出這個規則。一開始是在一個論壇上看到有imagej Macro的解法,幾行就可以生成jet lut

r=newArray(256);g=newArray(256);b=newArray(256); 
for (i=0;i<256;i++) { 
    i4=4*i/256; 
    r[i]=255*minOf(maxOf(minOf(i4-1.5,-i4+4.5),0),1); 
    g[i]=255*minOf(maxOf(minOf(i4-0.5,-i4+3.5),0),1); 
    b[i]=255*minOf(maxOf(minOf(i4+0.5,-i4+2.5),0),1); 

setLut(r,g,b); 


雖然可以生出jet Colormap,但是中間的i4-1.5或是-i4+4.5,弄不懂怎麼修改,所以暫時放下,繼續找其他文章。後來在stackoverflow的討論串裡找到解法。


這樣的解法就容易修改了,我只要縮窄原有的範圍(從256縮到151),再做位移(從0往後位移到25)就可以了。根據這位提問者的需求,做出來的就會是這樣。



在imagej點選File/New/Text Window,把下面這段貼到視窗裡,然後把Text Window的Language設定為IJ1 Macro,再按下Run/Run就可以修改Lut了。

在前幾行的xrange、xoffset和NanFill的數值都可以依照需求修改。
==============================================

run("8-bit");

xrange = 126; //灰階值範圍,不需特製的話,使用256
xoffset = 25; //最低灰階值,不需特製的話使用0
NaNFill = 255; //沒有灰階的部份填充的顏色,白色為255,黑色為0

r=newArray(256);
g=newArray(256);
b=newArray(256);

for (i=0;i<xrange;i++) {
x= -1 + i/(xrange/2);
r0 = 1.5 - abs(2*x -1);
    g0 = 1.5 - abs(2*x);
    b0 = 1.5 - abs(2*x +1);
   
    //clamp
    r[i+xoffset]=255*minOf(maxOf(r0,0),1);
    g[i+xoffset]=255*minOf(maxOf(g0,0),1);
    b[i+xoffset]=255*minOf(maxOf(b0,0),1);
}

//填充灰階值下界
for (i=0;i<xoffset;i++) {
r[i]=NaNFill;
    g[i]=NaNFill;
    b[i]=NaNFill;
}
//填充灰階值上界
for (i=xrange+xoffset;i<256;i++) {
r[i]=NaNFill;
    g[i]=NaNFill;
    b[i]=NaNFill;
}

setLut(r,g,b);  //設定LUT
run("Show LUT");  //顯示LUT的list



2019年2月5日

使用Imagej做影像分割的流程-使用PlantCV的流程教學

農業研究上,常常需要分析大量植物的表型,如果用傳統的個別測量方式,費時不少。因此目前已經有開發了自動機器人,固定為大量植栽攝影,而攝影所得的影像資料也是透過程式來加速處理影像分割和分析。這次介紹的網站裡的程式-PLANTCV就是一個用在此目的的PYTHON程式。
https://plantcv.readthedocs.io/en/latest/


我閱讀了其中流程的教學文件,覺得很實用,雖然是PYTHON的程式,不過影像分割的流程其實也能應用在imagej。因此本文就是將其Python流程轉換成Imagej流程的紀錄,如果是用imagej來處理大量影像資料,例如幾百張圖,那麼也只要把這些步驟錄製在macro裡再做修改就可以大量處理了。

以下是一些步驟的個別說明:

轉換色彩空間
以往我使用imagej做影像分割時,都是先用原始影像轉8 bits灰階,再用threshold或是原始彩色影像用color threshold來取得目標物。不過其實可以先將原始影像轉成HSV或是Lab的影像。因為有些目標物在RGB的channel之中未必會明顯突出,但是在HSV或是LAB影像中反而會變得明顯,因此轉成這類型影像後,再針對其中比較明顯的channel來做影像分割。
imagej的指令
Image/Type/RGB Stack
Image/Type/HSB Stack
Image/Type/Lab Stack


針對不同部位使用不同的影像空間的channel
由於植物並非色彩均質的物體,有些地方比較綠、有些地方黃甚至紅,因此可以使用不同色彩空間的channel來分割不同部位,最後再把這些組合在一起就可以。例如在Lab的色彩空間中,a channel是偏紅色綠色的,而b channel是偏藍色黃色的,所以就可以從這些部位個別取出要的部位再組合。


設定閾值做二值化影像
這部份就是直接設定Thershold,選擇要的部位來做二值化影像
imagej的指令
Image/Adjust/Threshold
Image/Adjust/Color Threshold




去除雜點、去噪
在二值化的影像中,難免會有一些噪點,可以使用中值濾波器,或是其他濾波器將雜點去除
imagej的指令
Process/Noise/Despeckle...
Process/Filters/Median..




影像邏輯運算
通常使用or運算做聯集,例如前述的過程得到一個以葉子為主的二值化影像,和一個以莖為主的二值化影像,接下來就可以用這兩個影像做聯集,也就是or運算,將兩個影像合併成一組。
imagej的指令
Process/Image Calculator


影像相減與相加
影像相減subtract,常用於減去背景,例如兩張圖片一張為沒有植株的盆栽,另一個是有植株的,兩者影像相減可得到植株的影像(但通常還需要再處理才能使用)
影像相加add,例如產生兩張圖,一張是邊緣強化的(透過sobel運算),另一張是增強對比的(經過laplace濾波),把兩張影像相加,可以得到兩者的優點,邊緣清楚且對比明顯
imagej的指令
Process/Image Calculator
Process/Find edge
Process/Enhance Contrast


從遮罩mask取得原圖目標物的選取區
例如有一個已經經過多次處理得到的二值化影像,在目標物的部份已經有了mask,可以將此mask轉成selection,再於原圖中套用此選取區
imagej的指令
Edit/Selection/Creation Selection(使用於二值化的影像,可得到選取區)
Edit/Selection/Restore selection(將其他影像的選取區套用於此)


以下用PlantCV的範例來做imagej的示範,共有兩個,並且用影片來說明。
---

---

流程:
1.原圖轉HSB --> 取得S channel影像 --> 設定閾值做二值化影像 --> 中值濾波去雜點  -->得到s_blu圖
2.原圖轉Lab --> 取得b channel影像 --> 設定閾值做二值化影像 --> 得到b圖
3.將s_blu圖和b圖做or運算,取得聯集後的影像,得到mask圖
4.將mask圖轉出選取區,套用在原圖上(或加入ROI manager)


第二個示範是紅外線攝影下的植物,有兩個影像,其中一個是有植物的,另一個是沒有植物的

--
--
流程:
1.原圖減去背景 -->  設定閾值做二值化影像  --->得到圖a
2.原圖做對比增強,利用laplace filter ,再用原圖減去這張圖,或是直接用Process/Enhance Contrast  --> 得到圖b
3.原圖做Find Edge尋找邊緣  --> 中值濾波  -->圖片反相 Edit/Invert   --> 得到圖c
4.圖b和圖c做add運算--> 設定閾值做二值化影像  ---> 侵蝕  --> 得到圖d
5.圖a和圖d做 or 運算 -->得到遮罩mask,再轉選取區套用在原圖上




形狀樣式(patten)的量化分析-使用imagej的套件 PAT-GEOM

最近看到一篇論文,是發展了一套Imagej的套件PAT-GEOM,可以用在分析生物身上的斑紋、花紋或條紋(以下這些統稱為樣式,patten)。

論文網址:
https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.13131

我仔細翻閱實作後,覺得是很有用的工具,以往我們在分析這些樣式的時候,通常也只分析面積有多大、顏色是什麼,再不然就只能說這個和那個「看起來」像,或是「看起來」比較亂...。其實不是我們分析的時候不想量化,而是根本舊不知道怎麼量化這些東西。我很高興能看看到這套工具,而且論文和guide都寫得很清楚,不僅能學習工具怎麼使用,也知道測量的原理和分析之後可以拿來做什麼,同時也提供了分析的例子,讓使用者能弄清楚工具的使用。

這套工具能檢測量化的項目有七項,包括樣式的形狀、方向性、尺寸、對比、分佈、整體樣式的分佈、隨機性等。在作者的網站上可以下載這個套件,也有使用說明和範例可以下載
http://ianzwchan.com/my-research/pat-geom/

安裝順序
1.安裝imagej(或是fiji)
2.將下載的套件解壓縮後放在imagej的plugins資料夾中
3.前往以下網址,下載橢圓傅立葉分析的套件。下載後一樣放在plugins資料夾
http://imagejdocu.tudor.lu/doku.php?id=plugin:analysis:fourier_shape_analysis:start
4.啟動imagej
5.在imagej工具列的Plugins最底下可以找到PAT-GEOM的選項


以下是這套工具能夠分析的項目說明,以及分析的例子。
除了以下文字說明外,我也錄製影片來說明分析方式和例子
---

---




一、形狀(Marking shape)
分析方式:利用橢圓傅立葉分析(Elliptical Fourier analysis)將樣式的形狀做量化。先將要分析的樣式圈出後成封閉的曲線,然後用橢圓傅立葉分析去用大量的函數擬合那個封閉曲線。這些擬合的函數都會有量化的數字,這些數字越接近,代表樣式越相似
例子:
比較杜鵑鳥和其宿主的蛋殼表面花紋
比較兩族群生物的樣式平均形狀,例如長頸鹿
鑑定有特殊色彩斑塊的個體(例如鯨鯊)
比較螃蟹甲殼上斑塊和背景的差異


二、樣式的方向性(Marking shape directionality)
分析方式:各個樣式的長寬比例和方向
例子:
比較食蚜蠅和黃蜂的身體條紋
比較動物身上的條紋和其背景的條紋(例如斑馬)
測量在基因調控下或是選汰壓力下的蝶翅眼點形狀
測量老虎的條紋的差異


三、尺寸(Marking size)
分析方式:除分析樣式的面積外,也會計算形心(幾何中心,Centroid)到樣式邊緣的距離平方總和,再平均之後取平方根的值(averaged centroid size)
例子:
比較人造的獵物和其模擬的對象差異(例如帝王蝶的捕食者實驗)
比較兩個不同族群的生物斑點尺寸的差異(例如七星瓢蟲的斑點)
比較動物和其背景的形狀尺寸


四、對比(Pattern contrast)
分析方式:使用變異係數(Coefficient of variation),用標準差除平均值的值。當對比越明顯,變異係數就越大
例子:
測定比目魚身上的色彩斑塊是否符合隨機取樣的背景基質
比較某一動物的兩個不同區域差異(例如會變色的花枝)


五、分佈(Marking distribution)
分析方式:像素陣列(Pixel matrix)。把生物的整個表面的樣式做成二值化影像(樣式處為1、其餘為0),成為Pixel matrix。把不同個體體表的樣式都做成Pixel matrix之後,可以疊加在一起成為單張圖,即可觀察「這群生物合起來的樣子長什麼樣」
例子:
將一整個族群的生物的「平均樣式」視覺化,例如把一群螃蟹的背甲做成pixel matrix看「這群螃蟹背甲長什麼樣」,觀察不同族群之間的差異
設計擬真的獵物,例如設計海參的警戒色樣式,取得一群海參的平均像素陣列後,觀察海參的平均樣式長什麼樣。


六、整體樣式分佈的方向性(Marking distribution directionality)
分析:分析體表全部的斑塊的角度和排列方向
例子:
在選汰壓力下,某個族群生物的斑塊是否會發展成線性排列(例如某些魚類體表的斑塊排列,或是蝶翅上的斑塊、眼斑)
比較兩種生物是否有相似的身體形狀

七、樣式的隨機性(Pattern randomness)
分析:將樣式存成GIF,用檔案的尺寸來表示隨機性。由於GIF檔案壓縮的特性,越隨機的圖樣,檔案尺寸會越大,因此可以用檔案的尺寸來代表隨機性
比較同種生物的不同型態的樣式(例如彩虹䗉螺Umbonium vestiarium此種生物具有多種不同的斑紋)
測定擬態的品質(例如杜鵑鳥和其宿主的蛋殼、岸蟹和其棲息的背景)

2019年2月3日

imagej中,進行二值化影像處理的一些方法

在imagej處理影像時,這個二值化影像處理的流程,我覺得是蠻實用的,包括形態學影像處理,如侵蝕、膨脹、開運算、閉運算、分水嶺、孔洞填充、邊界提取等。


在分析「影像中感興趣的位置的面積有多大」這樣的問題時,常常會遇到區塊之間彼此相連,無法分割計算,或是影像中有雜點、空洞,這時候我們需要進行影像的預處理來解決問題。
首先須將影像轉換成二值化影像,讓影像只有黑或白兩種顏色。在這之前要先決定「哪些顏色成為黑色,哪些成為白色」,所以需要一個界線數值,超過這個就變成黑(或白)。這個界線的數值稱為閾值,這個步驟在【Process/Binary/Make Binary】中進行

接下來需要一些形態學的處理,包含侵蝕、膨脹、開運算、閉運算。

侵蝕【Process/Binary/Erode】是將區塊的邊界向內收縮,所以小雜點就會消失。

膨脹Process/Binary/Dilate】是使邊界向外膨脹,可以填補區塊內的空洞。

開運算【Process/Binary/Open】,是先侵蝕再膨脹的過程,可以消除小的雜點,或是斷開兩個有細線相連的目標物。

閉運算【Process/Binary/Close】,是先膨脹再侵蝕,可以用來填充目標物內的細小孔洞,連接兩個鄰近的目標物。

提取影像的邊緣【Process/Binary/Outline】,實際的運算其實是把原圖減去侵蝕後的影像,如此可以得到影像的邊緣。

【Process/Binary/Watershed】,用在把兩個兩團相連的東西切開來,比方說細胞分裂時,中間會有牽絲的部份,這個指定可以斷開連結

【Process/Binary/Fill Holes】,顧名思義就是把影像中的洞填起來,將二值化影像中被黑色包圍的白色轉成黑色。

行道樹的接觸生長和死樹的比較

路過一棵行道樹的時候突然發現一個有趣的現象,活樹和死樹對於鐵絲的纏繞,會有不同的反應。

死樹不用說,當然就是凹下去被破壞。而活樹就不同了,右邊的活樹已經看不到鐵絲了,並不是不見了,而是樹把它包起來了!此謂「接觸生長

「當樹木碰到異物,像是其他或自己的樹幹、繩索、欄杆等物體,在接觸位置周圍的形成層會快速進行細胞分裂,以包覆異物,這會造成樹木養分和水分的運輸受阻,而成為樹木力學結構的缺陷。」




2019年2月1日

以光碟與跳棋黏土製作視錯覺的玩具

這幾年我們學校在暑假都會辦理給小學生參加的營隊,是由任教資優班的自然、數學老師設計課程來讓學生參加的。前幾年我們都是給學生做顯微鏡,以及使用uHandy做顯微觀察,不過今年想來玩別的,就想到可以玩玩視覺錯覺的活動。

就用光碟、跳棋和黏土來做平台吧,需要找一個能塞進光碟孔且還能卡住的跳棋,而黏土就是用來固定跳棋和光碟。此外,也利用黏土來平衡光碟。





平衡後且重心低的光碟可以轉得非常順,至少也能轉動30秒
有了這個平台,接下來就可以讓學生用發下的型紙做勞作了,分別讓學生做了這些:
1. 混色:在紙上用色筆畫出不同顏色,轉動後觀察顏色變化
2. 牛頓色盤:有了混色經驗,嘗試用不同顏色組合出白色
3. 錯覺旋轉盤:http://a-chien.blogspot.com/2009/08/blog-post_22.html
4. Benham Disk:http://faculty.washington.edu/chudler/benham.html









用氯化亞鈷試紙貼葉子的探究實驗

用氯化亞鈷試紙貼在葉子上,看上下表皮哪個先變色,這個其實是考題裡常常出現的題目,不過學生不一定有機會實作。其實是因為試紙通常在室內乾燥,拿到室外貼葉子恐怕都已經變紅色了。

不過因為我有了噹噹噹的電熱杯墊!!所以就可以帶學生出去外面做了,不過前提是室外也要有插座啦,叫學生帶著剪刀膠帶,我帶著電熱杯墊和氯化亞鈷試紙出去。

一人一片試紙,兩個學生一組,一人負責上表皮另一人負責下表皮,貼在同一片葉子的同樣地方,然後就是計時啦。

一種是計時三五分鐘後,看上下表皮變色的情形,另一種則是看下表皮多久會變色的時候,然後上表皮的狀況如何。

貼得好的話,上表皮幾乎是不會變色,除非是外界空氣滲入。

不過下表皮的變色通常有些問題,因為有浮起的葉脈,所以不容易用膠帶密貼,空氣會從邊緣滲入,所以變色都會從邊緣開始(請見下二圖)








後續還可以做什麼實驗呢?可以讓學生用摘下的葉子來做實驗,用POE的策略來進行課程。或是另外又把摘下的葉子插在水中進行實驗。透過這些實驗設計來看看實驗結果是否和預測的是否相同。
或者也可以測量葉面葉背的溫度,如果假設溫度和蒸散有關,那麼要怎麼排除因為葉面朝上、葉背朝下造成的影響。

又或者一些兩面都有氣孔的草本植物,以它們進行蒸散作用實驗,它們的葉子又會有怎麼樣的結果,溫度是否有會變化呢?

以手機進行定點攝影(3)-將地面觀測天空的照片和衛星圖結合

接續前篇,利用電腦下載手機定點攝影的照片後,可以做什麼事情呢?以往我都是拿來做縮時攝影,不過後來我想到一個應用。

曾經看過一個公民科學app的應用,那個app知悉你的位置之後,當有衛星通過你的上空拍下照片後,它會通知你,要你也往天空拍一張照片,然後就可以對照衛星圖和地面觀測的影像。

那麼我不如反過來做吧,就一直做地面觀測天空的攝影,然後再拿衛星圖來對照就好了啊。

先來看看成果。

這張組圖裡使用了中央氣象局的六種衛星觀測的產品,由左至右,由上而下分別是
真實色影像                   紅外線雲圖-黑白   紅外線雲圖-彩色
紅外線雲圖-色調強化  可見光雲圖            雷達合成回波圖
最右邊則是地面觀測的照片


這樣組合可以觀察什麼資訊呢?

  1. 真實色:可以觀察珊瑚紅色表示的低雲(或霧區)、灰、白色系為雲頂較高的雲層(詳細說明
  2. 色調強化:可以觀察雲頂的高度,藉此觀察對流系統的發展。顏色越綠等,代表雲頂溫度越低,也就是雲頂越高(詳細說明
  3. 雷達回波:有顏色反應的訊號代表降水的強度,藉此可以知道哪邊在降雨(詳細說明
  4. 紅外線黑白和可見光(直行第二排):兩張圖可以互相對照。可見光圖層可以看「雲的厚薄」,雲越厚就會越白。紅外線看的是「雲頂高度」,雲越高就會越冷,顏色就越白。
    在可見光看的到雲,但紅外線看不到雲,就是有厚厚的雲,但是雲的高度很低,例如冬天北部的層狀雲。
    在紅外線看得到雲,但可見光看不到雲,代表雲薄薄的,但是雲很高,那麼就是高空的卷雲、卷積雲等。(詳細說明

搭配地面觀測的照片,就可以從原本只是平面的地面觀測照片獲得多維度的資訊,同時看到這片雲的高度、厚度或是降水。我上面組圖特別選這張,是因為地面觀測到一片橫過天空的雲,而在衛星觀測上就看到一段橫過新竹到宜蘭的雲,剛好能互相對照。



以下這支程式就是用FTP下載手機照片後,再上中央氣象局網站抓取衛星圖,做成組圖的程式。預設都是晚上天黑後我才會執行,這樣才能夠下載到當天白天的所有圖片。
不過目前有幾個部份我還沒修正,
1.每天10:40會停止觀測,所以會少那時間的衛星資料,但是我程式沒有去處理這個部份,所以執行到這個時間的衛星圖組圖時,就會跳錯誤
2.有些手機照片會有拍攝錯誤的狀況,如果遇到這種照片,也會跳錯誤。

3.有兩個地方要手動改時間,第一個是下載衛星圖的區塊,如這段是下載2019年1月4日上午4點開始的衛星圖,程式中設定了下載此天從4點到20點的衛星圖。所以如果手機拍攝了兩天才下載,手機照片就會有兩天的圖片,而衛星圖就要下載兩天份的,所以就要改下載日期再重新執行一次。(不過氣象局只能下載48小時內的衛星圖)

start = datetime.datetime(2019, 1, 4, 4, 0, 0)  #下載日期

第二個要改的地方是程式末倒數十多行的地方,這個是用來製作組圖的。第一行是從上午六點開始抓存下的圖片來製作,不過中間會有幾個地方出錯。一是上午10:40的時候會錯,因為沒照片,再來就是有些衛星圖如果有錯無法讀取,也會在這時報錯。不過我一時也懶得改了

date = datetime.datetime(2019, 1, 4, 6, 0, 0)

以下是程式檔,我也同時放一份ipython notebook的在github


import subprocess
import datetime
import os
path = "/home/shawn-pc/桌面/SkyAndsatellite/"
if not os.path.isdir(path):  #檢查資料夾是否存在
    os.mkdir(path)

#從手機下載照片
from ftplib import FTP
HOST = "192.168.1.116"
PORT = 2221
pathFtp = "/opencamera"

ftp = FTP()
ftp.connect(HOST, PORT)
ftp.login()
ftp.cwd(pathFtp)
filenames = ftp.nlst()
for filename in filenames:
    if int(filename[13:15]) >=5 and int(filename[13:15]) <= 19:  #下載在5點到20點的照片
        dateFolder = filename[4:12]
        if not os.path.isdir(path + dateFolder+"/"):  #檢查日期資料夾是否存在
            os.mkdir(path + dateFolder+"/")

        if not os.path.isdir(path + dateFolder+"/"+"sky/"):  #檢查日期資料夾裡的sky資料夾是否存在
            os.mkdir(path + dateFolder+"/"+"sky/")


        local_filename = os.path.join(path + dateFolder+"/"+"sky/", filename)
        file = open(local_filename, 'wb')
        ftp.retrbinary('RETR '+ filename, file.write)
        print("download "+filename)
    
    file.close()
    ftp.delete(filename)

ftp.quit()


#下載衛星圖,先建立日期資料夾

import time

start = datetime.datetime(2019, 1, 4, 4, 0, 0)  #下載日期
pathDay = start.strftime(path+"%Y%m%d/")

if not os.path.isdir(pathDay):
    os.mkdir(pathDay)
print(pathDay)



#從cwb下載衛星圖
import requests

imgurl = "https://www.cwb.gov.tw/V7/observe/satellite/Data/{}/{}-{}.jpg"

def getImage(mapType):
    date = start
    
    pathSateType = pathDay + mapType +"/"
    if not os.path.isdir(pathSateType):
        os.mkdir(pathSateType)
    
    for i in range(6*(20-5)):
        mapTime = date.strftime("%Y-%m-%d-%H-%M")
        fname = mapType+"-"+mapTime+".jpg"
        with open(pathSateType+fname,'wb') as f:
            res = requests.get(imgurl.format(mapType,mapType,mapTime))
            f.write(res.content)
        date = date + datetime.timedelta(minutes = 10)

getImage("sbo") #可見光
getImage("s3p")  #紅外線雲圖彩色
getImage("s3o")  #紅外線雲圖黑白
getImage("s3q")  #色調強化
getImage("ts3p") #真實影像色彩



##下載雷達回波圖

radarUrl = "https://www.cwb.gov.tw/V7/observe/radar/Data/HD_Radar/CV1_TW_3600_{}.png"
def getRadar(date):
    pathRadar = pathDay + "radar/"
    if not os.path.isdir(pathRadar):
        os.mkdir(pathRadar)
        
    for i in range(6*(20-5)):
        mapTime = date.strftime("%Y%m%d%H%M")
        fname = "radar" + "-" + date.strftime("%Y-%m-%d-%H-%M") + ".png"
        
        with open(pathRadar+fname,'wb') as f:
            res = requests.get(radarUrl.format(mapTime))
            f.write(res.content)
        date = date + datetime.timedelta(minutes = 10)        


getRadar(date) #下載雷達回波



#製作組圖

from PIL import Image, ImageDraw, ImageFont
import re



"""檢查照片檔案是否存在"""
def is_jpg(filename):
    try:
        i=Image.open(filename)
        return i.format =='JPEG'
    except IOError:
        return False
    except FileNotFoundError:
        return False


"""取得衛星照片的路徑"""
def getMapPath(mapType,date):
    fname = mapType + date.strftime("-%Y-%m-%d-%H-%M.jpg")
    path = pathDay + mapType + "/" + fname
    return path


def getRadarPath(date):
    fname = "radar" + date.strftime("-%Y-%m-%d-%H-%M.png")
    path = pathDay + "radar" + "/" + fname    
    return path




"""取得天空拍照的照片路徑"""

def getSkyPath(date):
    imgNamePre =date.strftime("IMG_%Y%m%d_%H%M")
    for fname in os.listdir(pathDay + "sky/"):
        if re.match(imgNamePre+"\d+.jpeg", fname):
            return pathDay + "sky/" +fname

        


##製作組圖
def create_collage(width, height, listofimages, sky):
    cols = 4
    rows = 2
    thumbnail_width = width//cols
    thumbnail_height = height//rows
    size = thumbnail_width, thumbnail_height
    new_im = Image.new('RGB', (width, height))
    ims = []
    
    

    for p in listofimages:
        im = Image.open(p)
        im.thumbnail(size)
        ims.append(im)
    i = 0
    x = 0
    y = 0
    for col in range(cols):
        for row in range(rows):
            if i < len(ims):
                print(i, x, y)
                new_im.paste(ims[i], (x, y))
                i += 1
                y += thumbnail_height
        x += thumbnail_width
        y = 0

    #把天空攝影放在最右邊的欄
    
    imSky = Image.open(sky)
    area = (416,4,416+1060,4+3076)   #切割天空
    imSky = imSky.crop(area)    
    
    sizeSky = thumbnail_width , height
    imSky.thumbnail(sizeSky)

    x = thumbnail_width * 3

    new_im.paste(imSky, (x, 0))


    new_im.save(date.strftime(pathDay+"%Y%m%d%H%M.jpg"))
    #new_im.show()


date = datetime.datetime(2019, 1, 4, 6, 0, 0)
pathDay = date.strftime(path+"%Y%m%d/")
for i in range(90): s3o = getMapPath("s3o",date) s3p = getMapPath("s3p",date) s3q = getMapPath("s3q",date) sbo = getMapPath("sbo",date) ts3p = getMapPath("ts3p",date) radar = getRadarPath(date) sky = getSkyPath(date) #檢查路徑是否存在,如果存在才進行下一步 #how? listofimages=[ts3p,s3q,s3o,sbo,s3p,radar] create_collage(2000, 1000, listofimages,sky) date = date + datetime.timedelta(minutes = 10)

你可能對這有興趣

Related Posts Plugin for WordPress, Blogger...