data_folder = Path("/home/bathiany/Sebastian/datamining/edges/obsscan")
data_set = DataSet([data_folder / 'ice_conc_nh_ease2-250_cdr-v2p0_remapbilt200_September.nc'], 'ice_conc')

#dataset goes from 1979 - 2015

sigma_d = unit('120 km')   #200
sigma_t = unit('1 year')

gamma = 1e10
scaling_factor = gamma * unit('1 km/year')
sobel_delta_t = unit('1 year')
sobel_delta_d = sobel_delta_t * scaling_factor
sobel_weights = [sobel_delta_t, sobel_delta_d, sobel_delta_d]

"""
Explanation: Enter your settings here
End of explanation
"""


from datetime import date, timedelta

box = data_set.box

print("({:.6~P}, {:.6~P}, {:.6~P}) per pixel".format(*box.resolution))
for t in box.time[:3]:
    print(box.date(t), end=', ')
print(" ...")

dt = box.time[1:] - box.time[:-1]
print("time steps: max", dt.max(), "min", dt.min())

"""
Explanation: Load the data
End of explanation
"""


yearly_data_set = data_set
box = yearly_data_set.box
data = yearly_data_set.data

smooth_data = gaussian_filter(box, data, [sigma_t, sigma_d, sigma_d])

my_cmap = matplotlib.cm.get_cmap('rainbow')  #YlGnBu
my_cmap.set_under('w')

plot_orthographic_np(box, yearly_data_set.data[0], cmap=my_cmap, vmin=0.1) #1979


plot_orthographic_np(box, yearly_data_set.data[28], cmap=my_cmap, vmin=0.1) #2007

plot_orthographic_np(box, yearly_data_set.data[33], cmap=my_cmap, vmin=0.1) #2012


plot_orthographic_np(box, yearly_data_set.data[36], cmap=my_cmap, vmin=0.1) #2015


"""
Explanation: Smoothing
Smoothing in the lattitudinal direction (N-S) is not affected by the projection; distances do not vary since the circles of equal longitude are all greater circles on the sphere. We apply a Gaussian filter where we reflect data in the lattitudinal and time directions and wrap in the longitudinal direction. Next, to filter in the longitudinal direction, we need to use a different filter width for each lattitude.
End of explanation
"""


sb = sobel_filter(box, smooth_data, weight=sobel_weights)
pixel_sb = sobel_filter(box, smooth_data, physical=False)

"""
Explanation: Sobel filtering
The Sobel filter has the same problem as the Gaussian filter, but the solution is easier. We just correct for the magnitude of the Sobel response by multiplying the longitudinal component by the cosine of the latitude.
End of explanation
"""


signal = 1/sb[3]

## quartiles

signalarray = np.asarray(signal).reshape(-1)
signal_no0 = np.ma.masked_equal(signalarray,0)
signal_no0=signal_no0.compressed()

#np.percentile(signal_no0, [25, 50, 75, 100])

## There is no control data for observations => must choose thresholds based on data itself.
# Here, use 75 and 50%ile (quartiles 2 and 3)

upper_threshold=np.percentile(signal_no0, 95)
lower_threshold=np.percentile(signal_no0, 50)

from hyper_canny import cp_edge_thinning, cp_double_threshold

# use directions of pixel based sobel transform and magnitudes from calibrated physical sobel.
dat = pixel_sb.transpose([3,2,1,0]).astype('float32')
dat[:,:,:,3] = sb[3].transpose([2,1,0])
mask = cp_edge_thinning(dat)
thinned = mask.transpose([2, 1, 0])
dat = sb.transpose([3,2,1,0]).copy().astype('float32')

edges = cp_double_threshold(data=dat, mask=mask, a=1/upper_threshold, b=1/lower_threshold)
m = edges.transpose([2, 1, 0])

years = np.array([d.year for d in box.dates])
my_cmap = matplotlib.cm.get_cmap('rainbow')  #YlGnBu
my_cmap.set_under('w')
plot_orthographic_np(box, (years[:,None,None]*m).max(axis=0),  cmap=my_cmap, vmin=1979)

"""
Explanation: Determine hysteresis settings
End of explanation
"""


