import numpy as np
import xarray as xr
import sys
import os
import glob

import datetime

import argparse

def GetMaskedVariable(xd, name) :
    FillValue=1.e+20
    x=xd[name].data
    x[x>FillValue]=np.nan; x[x<-FillValue]=np.nan
    return x
#
#####################################################################
#
#
infile="etzh/qrain_etzh.nc"

indata=xr.open_dataset(infile)
print(indata)

print(indata["Times"].attrs["units"])
ff=indata["Times"].attrs["units"].replace("day as","")

date=[]
for t in indata["Times"].data :
    y=int(t/10000)
    m=int((t-y*10000)/100)
    d=int(t-(y*10000+m*100))
    date.append(datetime.datetime(y, m, d, 0, 0, 0)+datetime.timedelta(days=0.5))

print("Range of dates : ", date[0], date[-1])

lon=GetMaskedVariable(indata, "LON")
lat=GetMaskedVariable(indata, "LAT")
x=GetMaskedVariable(indata, "QRAIN")

outdata=xr.Dataset(data_vars=dict(
    Rainf=(["time","y","x"], x,
    {'units':indata["QRAIN"].attrs["units"], 'long_name':indata["QRAIN"].attrs["description"]})),
    coords=dict(
        lon=(["y","x"], lon,
             {'units':indata["LON"].attrs["units"], 'long_name':indata["LON"].attrs["description"]}),
        lat=(["y","x"], lat,
             {'units':indata["LAT"].attrs["units"], 'long_name':indata["LAT"].attrs["description"]}),
        time=("time", date)
        ),
    attrs=dict(description="NOAH-MP output from MeteoCat"),
    )

outdata.sel(time=slice('1989-01-01','2013-12-31')).to_netcdf("test_Rainf.nc")
outdata.close()