import numpy as np
import netCDF4 as nc
from pyhdf.SD import SD,SDC
from pyhdf.HDF import *
from pyhdf.VS import *
import glob
import h5py
import os

def hdfvsread(fname,varname):
  f = HDF(fname,SDC.READ)
  vs = f.vstart()
  vd = vs.attach(varname)
  temp = [data for data in vd[:]]
  test = [item[0] for item in temp]
  vd.detach()
  vs.end()
  f.close()
  return np.array(test)

def hdfSDread(fname,varname):
  f = SD(fname,SDC.READ)
  test = np.array(f.select(varname))
  f.end()
  return test

def ROIindex(Latitude,Longitude,Lat1,Lat2,Lon1,Lon2):
  LatM,Latm,LonM,Lonm = Lat1,Lat2,Lon1,Lon2
  if(Lat1<Lat2):
    LatM,Latm = Lat2,Lat1
  if(Lon1<Lon2):
    LonM,Lonm = Lon2,Lon1
  idx1,idx2 = np.where(Latitude > Latm)[0],np.where(Latitude < LatM)[0]
  idx3,idx4 = np.where(Longitude > Lonm)[0],np.where(Longitude < LonM)[0]
  select_index1,select_index2 = np.intersect1d(idx1,idx2),np.intersect1d(idx3,idx4)
  select_index = np.intersect1d(select_index1,select_index2)
  return select_index

def ROIConvectionindex(CloudType,ROI_idx):
  temp = set(np.where(CloudType == 8)[0])
  return np.intersect1d(np.array(sorted(temp)),ROI_idx)

def longest_consecutive_sequence(nums):
  if len(nums) == 0:
    return [],0
  nums = set(nums)
  longest_sequence = []
  current_sequence = []
  for num in nums:
    if num - 1 not in nums:
      current_sequence = [num]
      current_num = num + 1
      while current_num in nums:
        current_sequence.append(current_num)
        current_num += 1
      if len(current_sequence) > len(longest_sequence):
        longest_sequence = current_sequence
  return longest_sequence,len(longest_sequence)

def IsLand(latC,lonC):
  IsLand = 1 #Land
  if(latC - lonC > 19.0 and latC - lonC < 21.0):
    IsLand = 2 #Coastal Region
  elif(latC - lonC <= 19.0):
    IsLand = 0 #Ocean
  return IsLand

f = open('DeepConvective.txt','r')
lines = f.readlines()
f.close()
fname_list = []
for i in range(len(lines)):
  fname_list.append(lines[i][0:-1])

cldclass = '/aosc/eos20/syshan/CloudSat/cloudsat/cldclass/'
_2cice = '/aosc/eos20/syshan/CloudSat/cloudsat/2c_ice/'

fselect = []
for i in range(len(fname_list)):
  fname = fname_list[i]
  year = fname[0:4]
  #cldclassf = cldclass + str(year) + '/' + fname
  cldclassf = cldclass + fname
  print(cldclassf)

  Latitude,Longitude = hdfvsread(cldclassf,'Latitude'),hdfvsread(cldclassf,'Longitude')
  CloudLayer,CloudType = hdfvsread(cldclassf,'CloudLayer'),hdfSDread(cldclassf,'CloudLayerType')
  ROI_idx = ROIindex(Latitude,Longitude,-20,-40,-40,-60)
  LongestConvective,L = longest_consecutive_sequence(ROIConvectionindex(CloudType,ROI_idx))
  LatM,Latm = np.max(Latitude[LongestConvective]),np.min(Latitude[LongestConvective])
  if(LatM - Latm < 1.0 or LatM == np.max(Latitude[ROI_idx]) or Latm == np.min(Latitude[ROI_idx])):
    continue
  LonM,Lonm = np.max(Longitude[LongestConvective]),np.min(Longitude[LongestConvective])
  if(LonM == np.max(Longitude[ROI_idx]) or Lonm == np.min(Longitude[ROI_idx])):
    continue
  fselect.append(fname)

CTH_save,CBH_save,THICKNESS_save = [],[],[]
IWP_save,IWCHEIGHT_save = [],[]
LATC_save,LONC_save,DATE_save,IsLand_save = [],[],[],[]
for i in range(len(fselect)):
  fname = fselect[i]
  year,starttimeHH,starttimeMM,fdate = fname[0:4],int(fname[7:9]),int(fname[9:11]),fname[0:7]
  cldclassf = cldclass + fname
  Latitude,Longitude,Time = hdfvsread(cldclassf,'Latitude'),hdfvsread(cldclassf,'Longitude'),hdfvsread(cldclassf,'Profile_time')
  CloudLayer,CloudType = hdfvsread(cldclassf,'CloudLayer'),hdfSDread(cldclassf,'CloudLayerType')
  ROI_idx = ROIindex(Latitude,Longitude,-20,-40,-40,-60)
  LongestConvective,L = longest_consecutive_sequence(ROIConvectionindex(CloudType,ROI_idx))
  LatC,LonC,TimeC = np.mean(Latitude[LongestConvective]),np.mean(Longitude[LongestConvective]),round(np.mean(Time[LongestConvective])/60)
  print(LatC,LonC)
  LatM,Latm = LatC+0.5,LatC-0.5
  idx = np.intersect1d(np.where(Latitude[LongestConvective] > Latm)[0],np.where(Latitude[LongestConvective] < LatM)[0])
  select_idx = []
  for j in idx:
    select_idx.append(LongestConvective[int(j)])
  #print(select_idx)

  ##### Time #####
  if(TimeC > 60):
    HH,MM = int(TimeC/60),TimeC-60*int(TimeC/60)
  else:
    HH,MM = 0,TimeC
  timelocationMM,timelocationHH = MM+starttimeMM,HH+starttimeHH
  if(timelocationMM >= 120):
    timelocationMM,timelocationHH = timelocationMM-120,timelocationHH+2
  if(timelocationMM >= 60):
    timelocationMM,timelocationHH = timelocationMM-60,timelocationHH+1
  #print(fname[0:7],timelocationHH,timelocationMM)
  ##### Time #####

  ##### Cloud Thickness ######
  CTH,CBH = hdfSDread(cldclassf,'CloudLayerTop'),hdfSDread(cldclassf,'CloudLayerBase')
  cth,cbh = [],[]
  for j in select_idx:
    idx = np.where(CloudType[j,:] == 8.)[0]
    cth.append(CTH[j,idx][0])
    cbh.append(CBH[j,idx][0])
  cthmean,cbhmean = np.mean(cth),np.mean(cbh)
  thickness = cthmean - cbhmean
  CTH_save.append(cthmean) ## variable
  CBH_save.append(cbhmean) ## variable
  THICKNESS_save.append(thickness) ## variable
  #print(cthmean,cbhmean,thickness)
  ##### Cloud Thickness ######

  fid = fselect[i].split('_')[0] + '_' + fselect[i].split('_')[1]
  print(fid)

  ##### Cloud Ice #####
  fname = glob.glob(_2cice + fid + '*')[0]
  #print(fname)
  IWP,IWC,H = hdfvsread(fname,'ice_water_path'),hdfSDread(fname,'IWC'),hdfSDread(fname,'Height')
  iwp = np.nanmean(IWP[select_idx])
  IWP_save.append(iwp) ## variable
  iwcc = []
  for j in select_idx:
    iwcc.append(np.nansum(IWC[j] * H[j]) / np.nansum(IWC[j]))
  iwcheight = np.nanmean(iwcc)/1000.0
  IWCHEIGHT_save.append(iwcheight)
  ##### Cloud Ice #####

  ##### Ze #####
  #fname = glob.glob(geoprof + str(year) + '/' + fid + '*')[0]
  #print(fname)
  #####ze,H = hdfvsread(fname,'ze')
  ##### Ze #####

  LATC_save.append(LatC)
  LONC_save.append(LonC)
  DATE_save.append(fdate)
  IsLand_save.append(IsLand(LatC,LonC))

savename = './1DeepConvective.h5'
saveid = h5py.File(savename,'w')
dset = saveid.create_dataset('CTH',data=CTH_save,dtype=float,compression='gzip')
dset.attrs['Name'] = 'Cloud Top Height'
dset.attrs['Unit'] = 'km'
dset = saveid.create_dataset('CBH',data=CBH_save,dtype=float,compression='gzip')
dset.attrs['Name'] = 'Cloud Base Height'
dset.attrs['Unit'] = 'km'
dset = saveid.create_dataset('THICKNESS',data=THICKNESS_save,dtype=float,compression='gzip')
dset.attrs['Name'] = 'Cloud Thickness'
dset.attrs['Unit'] = 'km'
dset = saveid.create_dataset('IWP',data=IWP_save,dtype=float,compression='gzip')
dset.attrs['Name'] = 'Ice Water Path'
dset.attrs['Unit'] = 'g / m^2'
dset = saveid.create_dataset('IWCHEIGHT',data=IWCHEIGHT_save,dtype=float,compression='gzip')
dset.attrs['Name'] = 'Ice Water Content Centroid Height'
dset.attrs['Unit'] = 'km'
dset = saveid.create_dataset('LAT',data=LATC_save,dtype=float,compression='gzip')
dset.attrs['Name'] = 'Center Latitude'
dset = saveid.create_dataset('LON',data=LONC_save,dtype=float,compression='gzip')
dset.attrs['Name'] = 'Center Longitude'
dset = saveid.create_dataset('DAYNUM',data=DATE_save,dtype=int,compression='gzip')
dset.attrs['Name'] = 'Daynum'
dset = saveid.create_dataset('IsLand',data=IsLand_save,dtype=int,compression='gzip')
dset.attrs['Comment'] = '0:Ocean, 1:Land, 2:Coastal Region'
saveid.close()