接續
前篇，利用電腦下載手機定點攝影的照片後，可以做什麼事情呢？以往我都是拿來做縮時攝影，不過後來我想到一個應用。
曾經看過一個公民科學app的應用，那個app知悉你的位置之後，當有衛星通過你的上空拍下照片後，它會通知你，要你也往天空拍一張照片，然後就可以對照衛星圖和地面觀測的影像。
那麼我不如反過來做吧，就一直做地面觀測天空的攝影，然後再拿衛星圖來對照就好了啊。
先來看看成果。
這張組圖裡使用了中央氣象局的六種衛星觀測的產品，由左至右，由上而下分別是
真實色影像 紅外線雲圖-黑白 紅外線雲圖-彩色
紅外線雲圖-色調強化 可見光雲圖 雷達合成回波圖
最右邊則是地面觀測的照片
這樣組合可以觀察什麼資訊呢？
- 真實色：可以觀察珊瑚紅色表示的低雲(或霧區)、灰、白色系為雲頂較高的雲層（詳細說明）
- 色調強化：可以觀察雲頂的高度，藉此觀察對流系統的發展。顏色越綠等，代表雲頂溫度越低，也就是雲頂越高（詳細說明）
- 雷達回波：有顏色反應的訊號代表降水的強度，藉此可以知道哪邊在降雨（詳細說明）
- 紅外線黑白和可見光（直行第二排）：兩張圖可以互相對照。可見光圖層可以看「雲的厚薄」，雲越厚就會越白。紅外線看的是「雲頂高度」，雲越高就會越冷，顏色就越白。
在可見光看的到雲，但紅外線看不到雲，就是有厚厚的雲，但是雲的高度很低，例如冬天北部的層狀雲。
在紅外線看得到雲，但可見光看不到雲，代表雲薄薄的，但是雲很高，那麼就是高空的卷雲、卷積雲等。（詳細說明）
搭配地面觀測的照片，就可以從原本只是平面的地面觀測照片獲得多維度的資訊，同時看到這片雲的高度、厚度或是降水。我上面組圖特別選這張，是因為地面觀測到一片橫過天空的雲，而在衛星觀測上就看到一段橫過新竹到宜蘭的雲，剛好能互相對照。
以下這支程式就是用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)
