import warnings

import nc
import netCDF4
warnings.filterwarnings('ignore')
warnings.filterwarnings('ignore', category=DeprecationWarning)
import netCDF4
import pandas as pd
import numpy as np
from osgeo import gdal
import matplotlib.pyplot as plt
import math
import h5py
from pyhdf.SD import *
from osgeo import
 
 osr
import numpy as np

###
def array2raster(newRasterfn, rasterOrigin, xsize, ysize, array):
    """
     newRasterfn: 輸出tif路徑
     rasterOrigin: 原始柵格資料Extent
     xsize: x方向像元大小
     ysize: y方向像元大小
     array: 計算後的柵格資料
    """
    cols = array.shape[1]  # 矩陣列數
    rows = array.shape[0]  # 矩陣行數
    originX = rasterOrigin[0]  #
 
 起始像元經度
    originY = rasterOrigin[1]  # 起始像元緯度
    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Float32)
    # 括號中兩個0表示起始像元的行列號從(0,0)開始
    outRaster.SetGeoTransform((originX, xsize, 0, originY, 0, ysize))
    # 獲取資料集第一個波段，是從1開始，不是從0開始
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    outRasterSRS 
 
= osr.SpatialReference()
    # 程式碼4326表示WGS84座標
    outRasterSRS.ImportFromEPSG(4326)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()


###
from os import listdir
import re

inputf = 'E:/FSC_AVHRR/train_nc/'
output_TIF = 'E:/FSC_AVHRR/train_tif_4326/'
output_TIF_RE = 'E:/FSC_AVHRR/train_tif_6933/'
files = listdir(inputf)

pattern = '.nc$'
timer = 0

for i in files:
    temp_input = inputf + i

    temp_out_name = re.split(pattern, i)[0] + '.tif'
    temp_out_T_name = re.split(pattern, i)[0] + '_T.tif'
    temp_output = output_TIF + temp_out_name
    temp_output_T = output_TIF_RE + temp_out_T_name

    nc_data = netCDF4.Dataset(temp_input)
    FSC = nc_data.variables['scfg'][:]
    # plt.imshow(C)
    FSC_data = np.flipud(FSC[0, :, :])

    xsize = 0.05
    ysize = -0.05
    # 對於全球尺度rasterOrigin 為[-180(xmin),-90(ymin)]
    array2raster(temp_output, [-180, 90], xsize, ysize, FSC_data)
    option4 = gdal.WarpOptions(format='GTiff', srcSRS='EPSG:4326', dstSRS='EPSG:6933')
    FSC_RE = gdal.Warp(temp_output_T, temp_output, options=option4)
    timer = timer + 1
    print(timer)