Pythonを用いた空間統計分析のための衛星データ（LST, NDVI）の処理1

#!/usr/bin/env python
# coding: utf-8

import h5py
import pandas as pd
import tables
from pyhdf.SD import SD, SDC
import numpy as np
import pandas as pd
from scipy.interpolate import griddata
import seaborn as sns
import matplotlib.pyplot as plt
import os


##### 座標変換

### 座標変換後のCSVファイルを格納するフォルダを自動生成する。
new_dir_path_NDVI = './Data/NDVI_After_change_of_coordinate'
new_dir_path_LST = './Data/LST_After_change_of_coordinate'
os.makedirs(new_dir_path_NDVI, exist_ok=True)
os.makedirs(new_dir_path_LST, exist_ok=True)

### LSTの処理

# 座標変換前のh28v04データの取り出し

path = "./Data/LST_Before_change_of_coordinate/LST_h28v04/"
files = os.listdir(path)
ff = files[0:]

NDF = pd.read_csv(path + str(ff[0]),header=0,index_col=0)

for i in range(len(ff)):
    NDF = pd.read_csv(path + str(ff[i]),header=0,index_col=0)
    NDF00 = np.array(NDF)
    
    ### 座標投影処理
    # 例：h29v05の場合は北緯30-40、東経110-120である
    # h28v05はｈの1個ズレだから北緯30-40、東経100-110
    # h29v04はｖの1個ズレだから北緯40-50、東経110-120
    # Vは下側が小さい値（中心くらいが赤道で、そこから
    # 外れるほど値が大きくなるから）
    # https://gportal.jaxa.jp/gpr/assets/mng_upload/GCOM-C/
    # SGLI_Higher_Level_Product_Format_Description_jp.pdf
    # の3.7.2 タイル/EQA（全球EQA）参照
    C = np.ones((1,1200))
    #A = np.linspace(40,30,1200)
    A = np.linspace(50,40,1200)
    #B = np.linspace(110,120,1200)
    B = np.linspace(100,110,1200)
    A1 = np.reshape(A,(1200,1))
    AA = A1*C
    BB = 1/(np.cos(A1*np.pi/180))*B
    
    
    ### 欲しい座標部分の設定
    
    # 北海道全域
    ROI_Lat= [35.33 , 48.81]
    ROI_Lon= [135.53, 150.35]

    # パラメータ
    DDeg=10/1200

    CC2 = np.where((AA >= ROI_Lat[0]-0.05) & (AA <= ROI_Lat[1]+0.05)
     & (BB >= ROI_Lon[0]-0.05) & (BB <= ROI_Lon[1]+0.05))
    
    AAX = []
    BBY = []
    for g in range(len(CC2[0])):
        AA1 = AA[CC2[0][g]][CC2[1][g]]
        BB1 = BB[CC2[0][g]][CC2[1][g]]
        AAX.append(AA1)
        BBY.append(BB1)
        
    # 重複部分を消す
    AAAX = set(AAX)
    
    # mgrid(max,min,間隔)
    Latg = np.mgrid[max(AAAX):min(AAAX):-DDeg]
    Long = np.mgrid[min(BBY):max(BBY):DDeg]
    Latg,Long = np.meshgrid(Latg,Long)
    
    LST_XY = []
    for t in range(len(CC2[0])):
        #ND = ndvi[CC2[0][t]][CC2[1][t]]
        ND = NDF00[CC2[0][t]][CC2[1][t]]
        LST_XY.append(ND)
        
    #zi = griddata((BBY,AAX),Ndvi_XY,(Long,Latg),method='linear')
    zi = griddata((AAX,BBY),LST_XY,(Latg.T,Long.T),method='linear')
    
    Y = Latg.T
    X = Long.T
    # X,Yを1列にまとめる
    Y_qgis = Y.reshape(-1,)
    X_qgis = X.reshape(-1,)
    LST_qgis0 = zi.reshape(-1,)
    
    ck = np.c_[X_qgis,Y_qgis]
    dk = np.c_[ck,LST_qgis0]
    LST_qgis = pd.DataFrame(dk)
    LST_qgis = LST_qgis.rename({0:'X', 1:'Y', 2:'LST'},axis=1)
    # NDVIで欠損値のあるものを削除する
    LST_qgis = LST_qgis.dropna(how='any')
    LST_qgis.to_csv(new_dir_path_LST + "/LST_qgis_" + str(ff[i]))


################# NDVIの処理

# 座標変換前のh28v04データの取り出し

path = "./Data/NDVI_Before_change_of_coordinate/NDVI_h28v04/"
files = os.listdir(path)
ff = files[0:]

for i in range(len(ff)):
    NDF = pd.read_csv(path + str(ff[i]),header=0,index_col=0)
    NDF00 = np.array(NDF)
    
    ### 座標投影処理
    # 例：h29v05の場合は北緯30-40、東経110-120である
    # h28v05はｈの1個ズレだから北緯30-40、東経100-110
    # h29v04はｖの1個ズレだから北緯40-50、東経110-120
    # Vは下側が小さい値（中心くらいが赤道で、そこから外れるほど値が大きくなるから）
    # https://gportal.jaxa.jp/gpr/assets/mng_upload/GCOM-C/
    # SGLI_Higher_Level_Product_Format_Description_jp.pdf
    # の3.7.2 タイル/EQA（全球EQA）参照
    C = np.ones((1,1200))
    #A = np.linspace(40,30,1200)
    A = np.linspace(50,40,1200)
    #B = np.linspace(110,120,1200)
    B = np.linspace(100,110,1200)
    A1 = np.reshape(A,(1200,1))
    AA = A1*C
    BB = 1/(np.cos(A1*np.pi/180))*B
    
    
    ### 欲しい座標部分の設定
    
    # 関東地方のエリアを見たい時
    #ROI_Lat= [34.33 , 38.81]
    #ROI_Lon= [135.53, 141.35]

    # 全体を映したいとき
    ROI_Lat= [35.33 , 48.81]
    ROI_Lon= [135.53, 150.35]
    
    # 稚内エリア
    #ROI_Lat= [35.33 , 45.81]
    #ROI_Lon= [140.53, 147.35]

    # パラメータ
    DDeg=10/1200

    CC2 = np.where((AA >= ROI_Lat[0]-0.05) & (AA <= ROI_Lat[1]+0.05) 
    & (BB >= ROI_Lon[0]-0.05) & (BB <= ROI_Lon[1]+0.05))
    
    AAX = []
    BBY = []
    for g in range(len(CC2[0])):
        AA1 = AA[CC2[0][g]][CC2[1][g]]
        BB1 = BB[CC2[0][g]][CC2[1][g]]
        AAX.append(AA1)
        BBY.append(BB1)
        
    # 重複部分を消す
    AAAX = set(AAX)
    
    # mgrid(max,min,間隔)
    Latg = np.mgrid[max(AAAX):min(AAAX):-DDeg]
    Long = np.mgrid[min(BBY):max(BBY):DDeg]
    Latg,Long = np.meshgrid(Latg,Long)
    
    Ndvi_XY = []
    for t in range(len(CC2[0])):
        #ND = ndvi[CC2[0][t]][CC2[1][t]]
        ND = NDF00[CC2[0][t]][CC2[1][t]]
        Ndvi_XY.append(ND)
        
    #zi = griddata((BBY,AAX),Ndvi_XY,(Long,Latg),method='linear')
    zi = griddata((AAX,BBY),Ndvi_XY,(Latg.T,Long.T),method='linear')
    
    Y = Latg.T
    X = Long.T
    # X,Yを1列にまとめる
    Y_qgis = Y.reshape(-1,)
    X_qgis = X.reshape(-1,)
    NDVI_qgis0 = zi.reshape(-1,)
    
    ck = np.c_[X_qgis,Y_qgis]
    dk = np.c_[ck,NDVI_qgis0]
    NDVI_qgis = pd.DataFrame(dk)
    NDVI_qgis = NDVI_qgis.rename({0:'X', 1:'Y', 2:'NDVI'},axis=1)
    # NDVIで欠損値のあるものを削除する
    NDVI_qgis = NDVI_qgis.dropna(how='any')
    NDVI_qgis.to_csv(new_dir_path_NDVI + "/NDVI_qgis_" + str(ff[i]))