以手機進行定點攝影(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

https://github.com/ChihHsiangChien/CwbAndCellphone

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)