眼看着暑期實習即将結束,最近把近期的一些成果記錄一下,今後遇到類似的問題也可以進行參考。
開篇要說一下,這裡面需要用到osgeo的GDAL包,這個包可以從“https://www.lfd.uci.edu/~gohlke/pythonlibs/#pygame”下載下傳輪子進行安裝,我這裡用的是“GDAL-3.1.2-cp37-cp37m-win_amd64.whl”,3.1.2是版本,cp37表示适配python3.7的,大家可以根據自己的需求選擇下載下傳。
整個程式雖然是用于解析葵花8的AOD産品,但是其原理可以推廣至解析葵花8的其他産品乃至任意NetCDF格式資料。但是葵花8的AOD産品資料隻有三維(長*寬*變量),如果遇到更複雜的如同時包含海拔和時間維的資料就需要做一些改動了。
1.組成與功能介紹
程式由三部分組成,分别是“Download_AOD.py”、“AOD_NetCDF_to_GeoTIFF.py” 和 “Main.py”
“Download_AOD.py”主要用于下載下傳。下載下傳功能前一篇部落格已經寫過了,包括對AOD日均值資料和小時均值資料的下載下傳。當下載下傳日均值資料時,預設下載下傳昨天的AOD日均值資料。當下載下傳AOD小時均值資料時,預設下載下傳程式運作當天的中原標準時間9時至17時的逐小時均值資料,若程式運作時尚未更新完畢,則在将最新的資料下載下傳完成後自動終止運作。程式對正在下載下傳中的檔案以“.temp”格式在本地進行存儲,當目前檔案下載下傳完成後會轉為“.nc”格式存儲。在每次下載下傳開始前,程式會先将上次程式執行時未下載下傳完成的“.temp”檔案進行删除,并檢查存儲目錄中是否已經存在待下載下傳資料,若存在,則跳過,否則會開始下載下傳。
“AOD_NetCDF_to_GeoTIFF.py”用于将下載下傳後的AOD資料從NetCDF格式解析為TIFF格式,對于日均值和小時均值資料,解析字段名均為“AOT_L2_Mean”。根據官方文檔可知,NoData的栅格的值都設定為-32768,程式在解析時将這些值全部修改為0,也可以根據後期需求對NoData的值設定為其他特定的值。
“Main.py”用于運作。
2.存儲路徑與命名規則
假設程式運作時輸入的存儲路徑為footpath,則:
- 小時資料原始資料的存儲路徑為:footpath/YYYY-MM-DD/AOD_Hourly_Download
- 處理後的小時資料的存儲路徑為:footpath/YYYY-MM-DD/AOD_Hourly_Analysis
- 日資料原始資料的存儲路徑為:footpath/YYYY-MM-DD/AOD_Daily_Download
- 處理後的日資料的存儲路徑為:footpath/YYYY-MM-DD/AOD_Daily_Analysis
對于小時均值資料檔案,其原始檔案命名與官網保持一緻,格式為“H08_YYYYMMDD_hh00_1HARP030_FLDK.02401_02401.nc”,解析後的小時均值資料檔案命名格式為“YYYYMMDD_hh00.tif”。
對于日均值資料檔案,其原始檔案命名與官網保持一緻,格式為“H08_YYYYMMDD_0000_1DARP030_FLDK.02401_02401.nc”,解析後的日均值資料檔案命名格式為“YYYYMMDD_0000.tif”。
注:上述所有檔案夾與檔案命名中,YYYY表示四位數年份,MM表示兩位數月份,DD表示兩位數天,hh表示兩位數小時。檔案名中的時間均為UTC時間,使用時注意加8轉為中原標準時間。
3.代碼
3.1 下載下傳部分的代碼——“Download_AOD.py”
import os
import ftplib
import time
# 這個函數用于将日期從整型轉為FTP路徑所需的字元串
def getDateStr(yearNum, monNum, dayNum):
# 四位數年份
yearStr = str(yearNum)
# 兩位數月份
if monNum < 10:
monStr = "0" + str(monNum)
else:
monStr = str(monNum)
# 兩位數天
if dayNum < 10:
dayStr = "0" + str(dayNum)
else:
dayStr = str(dayNum)
return yearStr, monStr, dayStr
# 這個函數用于擷取前一天的日期号
def getYesterday(year, month, day):
if day == 1:
if month == 1:
year -= 1
month = 12
day = 31
elif month == 2 or month == 4 or month == 6 or month == 8 or month == 9 or month == 11:
month -= 1
day = 31
elif month == 5 or month == 7 or month == 10 or month == 12:
month -= 1
day = 30
elif month == 3:
# 閏年
if year % 4 == 0 and year % 400 == 0:
day = 29
month -= 1
# 閏年
elif year % 4 == 0 and year % 100 != 0:
day = 29
month -= 1
else:
day = 28
month -= 1
else:
day -= 1
return year, month, day
# 擷取檔案字尾名
def suffix(file, *suffixName):
array = map(file.endswith, suffixName)
if True in array:
return True
else:
return False
# 删除目錄下擴充名為.temp的檔案
def deleteFile(fileDir):
targetDir = fileDir
for file in os.listdir(targetDir):
targetFile = os.path.join(targetDir, file)
if suffix(file, '.temp'):
os.remove(targetFile)
class myFTP:
ftp = ftplib.FTP()
# 連接配接FTP,host是IP位址,port是端口,預設21
def __init__(self, host, port=21, YesdayNum=1):
self.ftp.connect(host, port)
self._dayNum = YesdayNum
# 登入FTP連接配接,user是使用者名,password是密碼
def Login(self, user, password):
self.ftp.login(user, password)
print(self.ftp.welcome) # 顯示登入資訊
# 下載下傳單個檔案,LocalFile表示本地存儲路徑和檔案名,RemoteFile是FTP路徑和檔案名
def DownLoadFile(self, LocalFile, RemoteFile):
bufSize = 102400
file_handler = open(LocalFile, 'wb')
print(file_handler)
# 接收伺服器上檔案并寫入本地檔案
self.ftp.retrbinary('RETR ' + RemoteFile, file_handler.write, bufSize)
self.ftp.set_debuglevel(0)
file_handler.close()
return True
# 下載下傳整個目錄下的檔案,LocalDir表示本地存儲路徑, emoteDir表示FTP路徑
def DownLoadFileTree(self, LocalDir, RemoteDir, choice):
# print("remoteDir:", RemoteDir)
# 如果本地不存在該路徑,則建立
if not os.path.exists(LocalDir):
os.makedirs(LocalDir)
# 擷取FTP路徑下的全部檔案名,以清單存儲
# 好像是亂序
self.ftp.cwd(RemoteDir)
RemoteNames = self.ftp.nlst()
RemoteNames.reverse()
# print("RemoteNames:", RemoteNames)
for file in RemoteNames:
# 先下載下傳為臨時檔案Local,下載下傳完成後再改名為nc4格式的檔案
# 這是為了防止上一次下載下傳中斷後,最後一個下載下傳的檔案未下載下傳完整,而再開始下載下傳時,程式會識别為已經下載下傳完成
Local = os.path.join(LocalDir, file[0:-3] + ".temp")
LocalNew = os.path.join(LocalDir, file)
'''
下載下傳小時檔案,隻下載下傳UTC時間1時至10時(中原標準時間9時至18時)的檔案
下載下傳的檔案必須是nc格式
若已經存在,則跳過下載下傳
'''
# 小時資料命名格式示例:H08_20200819_0700_1HARP030_FLDK.02401_02401.nc
if choice == 1:
if int(file[13:15]) >= 1 and int(file[13:15]) <= 10:
if not os.path.exists(LocalNew):
print("Downloading the file of %s" % file)
self.DownLoadFile(Local, file)
os.rename(Local, LocalNew)
print("The download of the file of %s has finished\n" % file)
elif os.path.exists(LocalNew):
print("The file of %s has already existed!\n" % file)
else:
pass
# 天資料命名格式示例:H08_20200802_0000_1DARP030_FLDK.02401_02401.nc
elif choice == 2:
if int(file[10:12]) == self._dayNum and not os.path.exists(LocalNew):
print("Downloading the file of %s" % file)
self.DownLoadFile(Local, file)
os.rename(Local, LocalNew)
print("The download of the file of %s has finished\n" % file)
elif int(file[10:12]) == self._dayNum and os.path.exists(LocalNew):
print("The file of %s has already existed!" % file)
self.ftp.cwd("..")
return
def close(self):
self.ftp.quit()
3.2 解析部分的代碼——“AOD_NetCDF_to_GeoTIFF.py”
# 子產品導入
import numpy as np
import netCDF4 as nc
from osgeo import gdal, osr, ogr
import os
import glob
def NC_to_tiffs(data, Output_folder):
# 讀取一下基本資訊
nc_data_obj = nc.Dataset(data)
Lon = nc_data_obj.variables["longitude"][:]
Lat = nc_data_obj.variables["latitude"][:]
# 讀取變量的時候,會自動根據scale factor對數值進行還原,但是Nodata的栅格會存儲為-32768
# 無論是日資料還是小時數居,變量名都是"AOT_L2_Mean"
AOD_arr = np.asarray(nc_data_obj.variables["AOT_L2_Mean"]) # 将AOD資料讀取為數組
# 這個循環将所有Nodata的值(即-32768)全部改為0
for i in range(len(AOD_arr)):
for j in range(len(AOD_arr[0])):
if AOD_arr[i][j] == -32768:
AOD_arr[i][j] = 0.0
# 影像的四秩
LonMin, LatMax, LonMax, LatMin = [Lon.min(), Lat.max(), Lon.max(), Lat.min()]
# 分辨率計算,其實可以寫死,都是2401*2401
N_Lat = len(Lat)
N_Lon = len(Lon)
Lon_Res = (LonMax - LonMin) / (float(N_Lon) - 1)
Lat_Res = (LatMax - LatMin) / (float(N_Lat) - 1)
# 此句代碼必須有
AOD_arr = np.array([AOD_arr])
for i in range(len(AOD_arr[:])):
# 建立.tif檔案
driver = gdal.GetDriverByName("GTiff")
TIFF_name = os.path.basename(data)
out_tif_name = Output_folder + "\\" + TIFF_name.split("_")[1] + "_" + TIFF_name.split("_")[2] + ".tif"
out_tif = driver.Create(out_tif_name, N_Lon, N_Lat, 1, gdal.GDT_Float32)
# 設定影像的顯示範圍
# -Lat_Res一定要是負的
geotransform = (LonMin, Lon_Res, 0, LatMax, 0, -Lat_Res)
out_tif.SetGeoTransform(geotransform)
# 擷取地理坐标系統資訊,用于選取需要的地理坐标系統
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326) # 定義輸出的坐标系為"WGS 84",AUTHORITY["EPSG","4326"]
out_tif.SetProjection(srs.ExportToWkt()) # 給建立圖層賦予投影資訊
# 資料寫出
out_tif.GetRasterBand(1).WriteArray(AOD_arr[i]) # 将資料寫入記憶體,此時沒有寫入硬碟
out_tif.FlushCache() # 将資料寫入硬碟
out_tif = None # 注意必須關閉tif檔案
3.3 執行程式——“Main.py”
import time
import os
import glob
import Download_AOD as Daod
import AOD_NetCDF_to_GeoTIFF as trans
dt = time.localtime(time.time())
_yearNum = dt.tm_year
_monNum = dt.tm_mon
_dayNum = dt.tm_mday
_yearStr = ""
_monStr = ""
_dayStr = ""
if __name__ == "__main__":
# 傳入IP位址
# 傳入的YesdayNum在下載下傳日資料時(昨天的)會需要
ftp = Daod.myFTP(host='ftp.ptree.jaxa.jp', YesdayNum=_dayNum - 1)
# 傳入使用者名和密碼,可以自行注冊并修改
ftp.Login('此處寫入您的賬号', '此處寫入您的密碼')
# 從目标路徑ftp_filePath将檔案下載下傳至本地路徑dst_filePath
dst_filePath = input("請輸入要存儲檔案的路徑:")
dst_filePath = dst_filePath + "/" + time.strftime("%F")
if not os.path.exists(dst_filePath):
os.makedirs(dst_filePath)
'''
下載下傳小時資料和日資料時,前置路徑都是:/pub/himawari/L3/ARP/030
下載下傳日資料時,示例路徑:/pub/himawari/L3/ARP/030/202008/daily/
下載下傳小時資料時,示例路徑:/pub/himawari/L3/ARP/030/202008/19/
'''
print("請選擇要下載下傳的資料:")
_choice = int(input("1.AOD小時資料(當天9:00至18:00) 2.AOD日均資料(昨天)\n"))
# Download_Path用于存儲下載下傳的原始資料
Download_Path = ""
# Analysis_Path用于存儲處理後的資料(即轉為TIFF後的資料)的檔案夾
Analysis_Path = ""
# 如果選擇為AOD小時資料
if _choice == 1:
_yearStr, _monStr, _dayStr = Daod.getDateStr(_yearNum, _monNum, _dayNum)
ftp_filePath = "/pub/himawari/L3/ARP/030" + "/" + _yearStr + _monStr + "/" + _dayStr + "/"
Download_Path = dst_filePath + "/AOD_Hourly_Download"
if not os.path.exists(Download_Path):
os.makedirs(Download_Path)
Daod.deleteFile(Download_Path) # 删除存儲路徑中的臨時檔案(也就是上次未下載下傳完整的檔案)
Analysis_Path = dst_filePath + "/AOD_Hourly_Analysis"
if not os.path.exists(Analysis_Path):
os.makedirs(Analysis_Path)
ftp.DownLoadFileTree(Download_Path, ftp_filePath, _choice)
# 如果選擇為AOD日資料(昨天的)
elif _choice == 2:
_yearNum, _monNum, _dayNum = Daod.getYesterday(_yearNum, _monNum, _dayNum)
_yearStr, _monStr, _dayStr = Daod.getDateStr(_yearNum, _monNum, _dayNum)
ftp_filePath = "/pub/himawari/L3/ARP/030" + "/" + _yearStr + _monStr + "/" + "daily" + "/"
Download_Path = dst_filePath + "/AOD_Daily_Download"
if not os.path.exists(Download_Path):
os.makedirs(Download_Path)
Daod.deleteFile(Download_Path) # 删除存儲路徑中的臨時檔案(也就是上次未下載下傳完整的檔案)
Analysis_Path = dst_filePath + "/AOD_Daily_Analysis"
if not os.path.exists(Analysis_Path):
os.makedirs(Analysis_Path)
ftp.DownLoadFileTree(Download_Path, ftp_filePath, _choice)
else:
print("選擇錯誤!")
# 下載下傳結束
ftp.close()
print("下載下傳完成!")
# 下面開始資料處理
# 讀取所有nc資料
data_list = glob.glob(Download_Path + "\\*.nc")
# for循環完成解析
for i in range(len(data_list)):
data = data_list[i]
trans.NC_to_tiffs(data, Analysis_Path)
print(data + "-----轉tif成功")
print("----轉換結束----")