使用Python自动完成Himawari-8(葵花8)卫星AOD数据下载与解析(转TIFF)

眼看着暑期实习即将结束,最近把近期的一些成果记录一下,今后遇到类似的问题也可以进行参考。

开篇要说一下,这里面需要用到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("----转换结束----")

 

 

相关推荐
©️2020 CSDN 皮肤主题: 游动-白 设计师:白松林 返回首页