2019年2月1日

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