In [13]:

    

import h5py


f0 = h5py.File('./rain/RAD_TF0005_A_20131104000000.h5')
p0 = f0.get('precipitation')[:]
p0[p0==-9999] = 0
f1 = h5py.File('./rain/RAD_TF0005_A_20131104010000.h5')
p1 = f1.get('precipitation')[:]
p1[p1==-9999] = 0
minx, maxx, miny, maxy = f1.attrs['grid_extent']
x = np.linspace(minx, maxx, num=p1.shape[1])
y = np.linspace(miny, maxy, num=p1.shape[0])


    

In [2]:

    

url = 'http://opendap.nationaleregenradar.nl/thredds/dodsC/radar/TF0100_A/2013/01/01/RAD_TF0100_A_20130101000000.h5'
import netCDF4


    

In [3]:

    

ds = netCDF4.Dataset(url)


    

In [4]:

    

x = ds.variables['east'][:]
y = ds.variables['north'][:]


    

In [5]:

    

P = ds.variables['precipitation']


    

In [6]:

    

p0 = P[:,:,7800].filled(0)
p1 = P[:,:,7801].filled(0)


    

In [18]:

    

fig, axes = plt.subplots(1,3, figsize=(13,4))
axes[0].imshow(p0, cmap='Blues', vmax=1, vmin=0)
axes[1].imshow(p1, cmap='Blues', vmax=1, vmin=0)
axes[2].imshow(p1 - p0, cmap='GnBu', vmax=0.1, vmin=-0.1)


    

    
Out[18]:





<matplotlib.image.AxesImage at 0x852cd90>






    

In [3]:

    

import cv2


    

In [108]:

    

N = lambda x:(x/x.sum())/(x/x.sum()).max()*255 
#N = lambda x:(x/x.max())*255
prvs = N(p0)
next = N(p1)
import urllib
import io
f = io.BytesIO(urllib.urlopen('https://upload.wikimedia.org/wikipedia/commons/thumb/5/55/Mona_Lisa_headcrop.jpg/375px-Mona_Lisa_headcrop.jpg').read())
p0 = matplotlib.image.imread(f, format='jpg').mean(2)
p1 = p0[:,::-1]
x = np.arange(p0.shape[1])
y = -np.arange(p0.shape[0])
flow = cv2.calcOpticalFlowFarneback(p0, p1, 0.7, 3, 11, 5, 5, 1.1, 0)
#flow = cv2.calcOpticalFlowFarneback(p0,next, 0.8, 5, 15, 10, 5, 1.2, 0)


    

In [110]:

    

fig, ax= plt.subplots()
ax.pcolorfast(x, y, p0, cmap='Blues')
ax.quiver(x[::10], y[::10], flow[::10,::10,0], flow[::10,::10,1], units='xy') #, linewidth=np.sqrt((flow[::10,::10]**2).sum(2))/5)


    

    
Out[110]:





<matplotlib.quiver.Quiver at 0x310280d0>






    

In [124]:

    

def warp_flow(img, flow, factor):
    h, w = flow.shape[:2]
    flow = -flow.copy()*factor
    flow[:,:,0] += np.arange(w)
    flow[:,:,1] += np.arange(h)[:,np.newaxis]
    res = cv2.remap(img, flow, None, cv2.INTER_LINEAR)
    return res


    

In [130]:

    

steps = np.linspace(0,10,num=6)
fig, axes = plt.subplots(1, steps.shape[0]+1, figsize=(steps.shape[0]*4, 3))
for i, step in enumerate(steps):
    axes[i].imshow(warp_flow(p0, flow, step),cmap='Blues')
    axes[i].set_title(step)
axes[-1].imshow(p1, cmap='Blues')


    

    
Out[130]:





<matplotlib.image.AxesImage at 0x3c9875d0>






    

In [45]:

    




    

In [ ]: