'''

comp = dict(zlib=True, complevel=1)

encoding = {var: comp for var in xds.data_vars}

xds.to_netcdf(fn, encoding=encoding)

result = np.zeros_like(arr, dtype=np.float32)

x=range(arr.shape[1])

y=range(arr.shape[2])

for i in x:

for j in y:

line_stack = pd.DataFrame(data=arr[:,i,j], dtype=np.float32,columns=['a'])

line_stack['a'].interpolate(method=method, axis=0, inplace=True, limit=limit)

if line_stack.a.isna().sum() >0:

line_stack['a']=-999

result[:, i, j] = line_stack['a'].values.astype(np.float32)

yco=spco[0]

xco=spco[1]

my_func_name='ApplyMannKendall'

print(datetime.now().strftime("%H:%M:%S")+' '+my_func_name)

xds_var=xds[var]

arr=np.array(xds_var)

#arr[arr<260]=np.nan

#pdb.set_trace()

#arr=nan_along_time_to_zero(arr) Doesnt work on Geofarm

#arr=interpolate_nan(arr,method='linear',limit=1)

arr=interpolate_along_time(arr, 'linear', 9)

result = np.apply_along_axis(mk.original_test, 0, arr) # Default alpha of 0.05

xds_res=xr.Dataset(coords={yco:xds.coords[yco],xco:xds.coords[xco]},

attrs=xds.attrs)

res_fields=['trend','h','p','z','Tau','s','var_s','slope','intercept']

i=0

for res in result:

field=res_fields[i]

if field=='trend':

res[res=='no trend']='0'

res[res=='increasing']='1'

res[res=='decreasing']='-1'

if field=='h':

res[res=='False']='0'

res[res=='True']='1'

res=res.astype(float)

xds_res[field]=([yco,xco],res)

i=i+1

'''

with fiona.open(path, 'r') as src:

for i, feature in enumerate(src):

if feature['properties']['id'] == ID:

index_i = i

return index_i

if not 1 <= day_of_year <= 365:

raise ValueError("Day of year must be between 1 and 365")

base_date = datetime(datetime.now().year, 1, 1)

target_date = base_date + timedelta(days=day_of_year - 1)

files_missing=[path_missing+f for f in os.listdir(path_missing)]

big_df=pd.DataFrame()

for file in files_missing:

big_df=big_df.append(pd.read_csv(file))

print('Mask out platform '+str(pf))

xds['pf_max']=((xds['platform'])/10000).astype(int)

xds=xds.where(xds.pf_max!=pf)

print('Correct daytime')

params=pd.read_csv(path_stats+study_area+'_daytime_params.csv')

param_a=float(params[params.month==month]['a'])

param_b=float(params[params.month==month]['b'])

corr_arr=dtc(t, param_a, param_b)

print('Calculating Decadal Anomalies')

'''

path_out= '/nfs/IGARSS_2022/Results_Laura/Results/New_Test/'

poly_lst = np.unique([os.listdir(path_out)[d].split('_')[0] for d in range(len(os.listdir(path_out)))])

path_missing='/nfs/IGARSS_2022/Results_Laura/to_process/'

path_stats='/nfs/IGARSS_2022/Results_Laura/Stats/'

#shp='D:/TIMELINE_SST/GIS/coast_dk/coast_dk.shp'

#shp = hard_drive + ':/TIMELINE_SST/GIS/sst_analysis_polygons/intersting_sst_analysis.shp'

shp = '/nfs/IGARSS_2022/Results_Laura/Shps/grid.shp'

#i = find_feature_index_by_id(shp, 3166)

start_date=date(2000,1,1)

end_date=date(2000,12,31)

aux=tl3_analysis_toolbox.aux()

tl_crop=tl3_analysis_toolbox.crop()

prep=tl3_analysis_toolbox.l3_lst_sst_prep()

reproj=tl3_analysis_toolbox.reproj()

shapefile=fiona.open(shp, "r")

'''

#missing_poly_lst=get_from_files(path_missing)

#missing_poly_lst=['3608']

#for i in range(len(shapefile)):

poly=[poly for poly in shapefile if poly['properties']['id']==int(poly_id)][0]

#poly_id= int(poly['properties']['id'])

'''

tiles = tl_crop.what_tiles_for_poly(poly)

print('Processing PolyID: ' + str(poly_id))

temp_res = 'Decades'

if temp_res == 'Daily':

date_range = pd.date_range(start_date, end_date, freq = 'd')

if temp_res == 'Decades':

daily_range = pd.date_range(start_date, end_date, freq = 'd')

date_range = daily_range[daily_range.day.isin([8,18,28])]

# make sure that each decade file is only processed once

for dt in date_range:

print('Calculating Anomalies for: ' + str(dt))

doi=dt.timetuple().tm_yday

xds_list=[]

if len(tiles)>0:

for tile in tiles:

#xds_mns=xr.Dataset(coords={'x':[],'y':[],'doi':[]})

fns = []

for year in range(1990,2023,1):

if temp_res == 'Decades':

if year > 2018:

version = 't01.40'

else:

version = 'v01.01'

date_f=str(year)+str(dt.month).zfill(2)+str(dt.day).zfill(2)

fn = aux.get_l3_file_from_ida(date_f, temp_res, 'SST', version, tile)

if fn !=None:

fns.append(fn)

else:

date_f=str(year)+str(dt.month).zfill(2)+str(dt.day).zfill(2)

fn=aux.get_l3_file_from_ida(date_f,temp_res,'SST','v01.01',tile)

if fn !=None:

fns.append(fn)

# Stack all Observations through the time axis

xds = xr.open_mfdataset(fns, combine='nested', concat_dim=['t'],

chunks={'x': 2300, 'y': 3250},

preprocess=prep.add_date)

xds=xds.drop(['lambert_azimuthal_equal_area'])

#if np.any(~np.isnan(xds['sst'][0].values)):

if xds.x.size > 0:

# Crop to Subpolygon Shape (1x1 Grad)

#xds=tl_crop.crop_ds_with_shp_tiles(xds,poly,crop,tile)

# Create binary mask from study area shape

coord_upper_left=(np.array(xds.x)[0],np.array(xds.y)[0])

resolution=1000

epsg_poly='4326'

epsg_mask='3035'

mask_shape=(xds.y.shape[0],xds.x.shape[0])

mask=tl_crop.create_mask_from_shp(poly,mask_shape,coord_upper_left,resolution,epsg_poly,epsg_mask)

xds=tl_crop.mask_ds_with_shp(xds, mask, poly, crop=True)

xds_list.append(xds) # append to list

#

chunksize=[1000,1000,37] # wird nicht mehr benötigt

if temp_res == 'daily':

var=['sst','view_time'] # keep only certain variables

# calculate mean and daily anomalies

xds_merged= reproj.mosaic_vars(xds_list,[var],chunksize)

xds_merged['sst_mean']=xds_merged['sst'].mean(dim='t')

xds_merged['sst_dif']=xds_merged['sst']-xds_merged['sst_mean']

xds_dif=xds_merged[['sst_dif']]

if temp_res == 'Decades':

'''

var=['qual_'+stats,

'sst_'+stats, 'sst_std',

'view_time_'+stats, 'valid_obs_count','platform']

if stats=='minmax':

var=['qual_min', 'qual_max',

'sst_min', 'sst_std', 'sst_max',

'view_time_min', 'view_time_max','valid_obs_count','platform']

xds_merged= reproj.mosaic_vars(xds_list,[var],chunksize)

# Apply daytime correction

if stats=='minmax':

stats1='max'

cor_arr=get_daytime_correction(np.array(xds_merged['view_time_'+stats1]),study_area,dt.month,path_stats)

xds_merged['sst_'+stats1]=xds_merged['sst_'+stats1]-cor_arr

stats2='min'

cor_arr=get_daytime_correction(np.array(xds_merged['view_time_'+stats2]),study_area,dt.month,path_stats)

xds_merged['sst_'+stats2]=xds_merged['sst_'+stats2]-cor_arr

# SST minmax = average between sst_min and sst_max

xds_merged['sst_'+stats]=(xds_merged.sst_min+xds_merged.sst_max)/2

else:

cor_arr=get_daytime_correction(np.array(xds_merged['view_time_'+stats]),study_area,dt.month,path_stats)

xds_merged['sst_'+stats]=xds_merged['sst_'+stats]-cor_arr

# Mask out NOAA-11

#xds_merged=mask_out_pf(xds_merged, 11)

# Mask out early and late observation times

#xds_merged=xds_merged.where((xds_merged['view_time_'+stats]>10)&(xds_merged['view_time_'+stats]<15))

#xds_merged['sst_mean_dec']=xds_merged['sst_mean'].mean(dim='t')

#xds_merged['sst_max_dec']=xds_merged['sst_max'].mean(dim='t')

#xds_merged['sst_min_dec']=xds_merged['sst_min'].mean(dim='t')

xds_merged['sst_med_dec']=xds_merged['sst_'+stats].mean(dim='t')

#xds_merged['sst_dif_mean']=xds_merged['sst_mean']-xds_merged['sst_mean_dec']

#xds_merged['sst_dif_min']=xds_merged['sst_min']-xds_merged['sst_min_dec']

#xds_merged['sst_dif_max']=xds_merged['sst_max']-xds_merged['sst_max_dec']

xds_merged['sst_dif_med']=xds_merged['sst_'+stats]-xds_merged['sst_med_dec']

xds_dif=xds_merged[['sst_dif_med']]

#xds_dif=xds_merged[['sst_dif_max']]

# Mask values with unrealistic anomalies

#xds_dif = xds_dif.where(xds_dif['sst_dif_max'] < 4)

#xds_dif = xds_dif.where(xds_dif['sst_dif_max'] > -4)

xds_dif = xds_dif.where(xds_dif['sst_dif_med'] < 4)

xds_dif = xds_dif.where(xds_dif['sst_dif_med'] > -4)

# Obs Count

#xds_dif['obs_count']=obs_count

xds_dif['obs_count'] = xds_merged['valid_obs_count']

#outfile=path_out+str(poly_id)+'_'+str(doi).zfill(3)+'_' + temp_res +'_dif.nc'

print('Load NC')

xds_dif=xds_dif.load()

print('Write Dif Netcdf')

outfile=path_out+str(poly_id)+'_'+str(doi).zfill(3)+'_' + temp_res +'_dif_'+stats+'.nc'

write_nc(xds_dif, outfile)

#xds_doi=xr.Dataset(coords={'x':xds_merged.coords['x'],'y':xds_merged.coords['y'], 'doi':doi})

xds_dif.close()

xds_merged.close()

xds.close()

gc.collect() # collect garbage

else:

print('Calculating monthly anomalies')

'''

path_in = '/nfs/IGARSS_2022/Results_Laura/Results/New_Test/'

path_out= '/nfs/IGARSS_2022/Results_Laura/Monthly_Results/New_Test/'

dec_doi = np.unique([os.listdir(path_in)[d].split('_')[1] for d in range(len(os.listdir(path_in)))])

poly_lst = np.unique([os.listdir(path_in)[d].split('_')[0] for d in range(len(os.listdir(path_in)))])

poly_lst_out = np.unique([os.listdir(path_out)[d].split('_')[0] for d in range(len(os.listdir(path_out)))])

poly_lst_out=[]

'''

if str(poly_id) in poly_lst_out:

print('Monthly Anomalies for PolyID ' + str(poly_id) + ' already calculated')

else:

print('Processing PolyID: ' + str(poly_id))

doi_perMonth = [dec_doi[i:i+3] for i in range(0, len(dec_doi), 3)]

for m in range(12):

#fns = [path_in + poly_id + '_' + str(doi_perMonth[m][i]).zfill(3) + '_Decades_dif.nc' for i in range(3)]

fns = [path_in + str(poly_id) + '_' + str(doi_perMonth[m][i]).zfill(3) + '_Decades_dif_'+stats+'.nc' for i in range(3)]

# 0 durch m ersetzen!

print(fns)

# Stack all Observations through the time axis

xds = xr.open_mfdataset(fns, combine='nested', concat_dim=['t'],

chunks={'x': 2300, 'y': 3250})

#xds_max = xds[['sst_dif_max']]

xds_max = xds[['sst_dif_med']]

xds_m = xds_max.groupby('t.year').mean('t')

xds_m['obs_count'] = xds['obs_count'].groupby('t.year').sum()

#xds_m = xds_m.rio.write_crs(3035)

print('Write Dif Netcdf')

outfile=path_out+str(poly_id)+'_'+str(m+1).zfill(2)+'_monthly_dif_'+stats+'.nc'

write_nc(xds_m, outfile)

# Apply MannKenddal

spco=['y','x']

#xds_res=ApplyMannKendall(xds_m,'sst_dif_max',spco)

xds_res=ApplyMannKendall(xds_m,'sst_dif_med',spco)

#xds_res = xds_res.rio.write_crs(3035) # set spatial reference

outfile=path_out+str(poly_id)+'_'+str(m+1).zfill(2)+'_monthly_mk_'+stats+'.nc'

write_nc(xds_res, outfile)

xds.close()

xds_res.close()

xds_m.close()

'''

parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)

#parser.add_argument('--start', help='No help',required=False)

#parser.add_argument('--stop', help='No help',required=False)

parser.add_argument('--stats', help='No help',required=False)

parser.add_argument('--study_area', help='No help',required=False)