In [1]:
import wradlib
import numpy as np
import os
import datetime as dt
In [2]:
%matplotlib
In [3]:
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib.collections import PatchCollection
import datetime
import warnings
warnings.simplefilter('once', DeprecationWarning)
In [4]:
from scipy import ndimage as ndi
from skimage import feature
from skimage.feature import match_template
In [5]:
import h5py
import pandas as pd
In [30]:
import matplotlib
from matplotlib.patches import Circle, Wedge, Polygon, Rectangle
from matplotlib.collections import PatchCollection
In [7]:
from skimage import measure
from skimage import filters
from scipy import ndimage
from skimage.measure import label, regionprops
import math
from matplotlib.patches import Ellipse
In [8]:
for_year = "2016"
rootdir = r"e:\data\radolan\ry"
tmpdir = r"e:\data\radolan\tmp"
h5file = 'ry_%s.hdf5' % for_year
hourlyfile = 'hdf/ry_hourly_%s.hdf5' % for_year
hourlyobjfile = 'hdf/ry_hourly_objects_%s.pickle' % for_year
tstart = "%s-01-01" % for_year
tend = "%s-12-31" % for_year
dffile = "exc_%s.csv" % for_year
nx = 900
ny = 900
thresh = 20. # mm/h
minarea = 10
maxarea = 1500
In [9]:
days = wradlib.util.from_to(dt.datetime.strptime(tstart, "%Y-%m-%d"),
dt.datetime.strptime(tend, "%Y-%m-%d"), tdelta=3600*24)
dtimes = wradlib.util.from_to(days[0].strftime("%Y-%m-%d 00:00:00"),
(days[-1]+dt.timedelta(days=1)).strftime("%Y-%m-%d 00:00:00"), tdelta=60*60)
hrs = np.arange(24).astype("i4")
In [10]:
dummy = regionprops(np.ones((4,4)).astype("i4"), intensity_image=np.ones((4,4)) )
keys = list(dummy[0])
props = [ dummy[0].__getitem__(key) for key in keys ]
keys.insert(0, "dtime")
props.insert(0, "1900-01-01 00:00:00")
In [91]:
df = pd.DataFrame( dict([(key, [props[i]]) for i,key in enumerate(keys)] ) )
In [92]:
with h5py.File(hourlyfile, 'r') as f:
for day in days:
print(day.strftime("%Y/%m/%d"), end="")
try:
dset = f[day.strftime("%Y/%m/%d")][:]
except KeyError:
print(" does not exist.")
continue
found = 0
for i, hr in enumerate(hrs):
hset = dset[i]
label_im = measure.label(hset > thresh, background=0)
nb_labels = len(np.unique(label_im))
regions = regionprops(label_im, intensity_image=hset)
for region in regions:
if (region.area < minarea) or (region.area > maxarea):
continue
found += 1
thetime = day.strftime("%Y-%m-%d") + " %02d:00:00" % hr
theprops = [region.__getitem__(prop) for prop in region]
theprops.insert(0, thetime)
df = df.append(dict([(key, theprops[i]) for i,key in enumerate(keys)] ), ignore_index=True)
print(" Found %d regions." % found)
In [11]:
#df.to_pickle(hourlyobjfile)
df = pd.read_pickle(hourlyobjfile)
In [32]:
df = df.set_index(pd.DatetimeIndex(df['dtime']))
In [13]:
df.keys()
Out[13]:
In [118]:
expandby = np.arange(50)
toleft = np.zeros((len(df), len(expandby))) * np.nan
toright = toleft.copy()
tobottom = toleft.copy()
totop = toleft.copy()
In [119]:
with h5py.File(hourlyfile, 'r') as f:
for i in range(len(df)):
dtime = dt.datetime.strptime(df.dtime.iloc[i], "%Y-%m-%d %H:%M:%S")
print(dtime)
try:
hset = f[dtime.strftime("%Y/%m/%d")][dtime.hour]
except KeyError:
continue
left, bottom, right, top = df.bbox.iloc[i][0], df.bbox.iloc[i][1], df.bbox.iloc[i][2], df.bbox.iloc[i][3]
for j, step in enumerate(expandby):
try:
toleft[i,j] = np.nanmean(hset[(left-step):right, bottom:top])
except IndexError:
continue
try:
toright[i,j] = np.nanmean(hset[left:(right+step), bottom:top])
except IndexError:
continue
try:
tobottom[i,j] = np.nanmean(hset[left:right, (bottom-step):top])
except IndexError:
continue
try:
totop[i,j] = np.nanmean(hset[left:right, bottom:(top+step)])
except IndexError:
continue
In [134]:
leftnorm = toleft / toleft[:,0].reshape((-1,1))
leftnorm[leftnorm>1] = np.nan
rightnorm = toright / toright[:,0].reshape((-1,1))
rightnorm[rightnorm>1] = np.nan
bottomnorm = tobottom / tobottom[:,0].reshape((-1,1))
bottomnorm[bottomnorm>1] = np.nan
topnorm = totop / totop[:,0].reshape((-1,1))
topnorm[topnorm>1] = np.nan
In [149]:
print("left")
for i, item in enumerate(leftnorm):
plt.plot(expandby, np.ma.masked_invalid(item), "b-", alpha=0.005)
print("right")
for i, item in enumerate(rightnorm):
plt.plot(expandby, np.ma.masked_invalid(item), "r-", alpha=0.005)
print("bottom")
for i, item in enumerate(bottomnorm):
plt.plot(expandby, np.ma.masked_invalid(item), "g-", alpha=0.005)
print("top")
for i, item in enumerate(topnorm):
plt.plot(expandby, np.ma.masked_invalid(item), "k-", alpha=0.005)
In [72]:
# Thresholds of hourly rainfall depths for detection of contiguous regions
threshs = np.arange(19,0,-1)
In [73]:
def get_regions(im, thresh):
"""Extract regions from im which exceed thresh.
"""
label_im = measure.label(im > thresh, background=0)
nb_labels = len(np.unique(label_im))
regions = regionprops(label_im, intensity_image=im)
return(regions)
In [74]:
#means = np.load("hdf/means_2016.numpy.npy")
#areas = np.load("hdf/areas_2016.numpy.npy")
In [75]:
means = np.zeros( (len(df), len(threshs)+1) )
areas = np.zeros( (len(df), len(threshs)+1) )
_dtime = None
with h5py.File(hourlyfile, 'r') as f:
for i in range(len(df)):
dtime = dt.datetime.strptime(df.dtime.iloc[i], "%Y-%m-%d %H:%M:%S")
bottom, left, top, right = df.bbox.iloc[i][0], df.bbox.iloc[i][1], df.bbox.iloc[i][2], df.bbox.iloc[i][3]
means[i, 0] = df.mean_intensity.iloc[i]
areas[i, 0] = df.area.iloc[i]
if dtime != _dtime:
print("")
print(dtime, end="")
# Only process new hourly set for a new datetime
try:
hset = f[dtime.strftime("%Y/%m/%d")][dtime.hour]
except KeyError:
continue
threshregions = [get_regions(hset, thresh) for thresh in threshs]
else:
print(".", end="")
for trix, tr in enumerate(threshregions):
for r in tr:
# Looking for region that contains core region
if (left >= r.bbox[1]) and \
(right <= r.bbox[3]) and \
(bottom >= r.bbox[0]) and \
(top <= r.bbox[2]):
# Found
means[i,trix+1] = r.mean_intensity
areas[i,trix+1] = r.area
_dtime = dtime
In [76]:
np.save("hdf/means_2016", means)
np.save("hdf/areas_2016", areas)
In [121]:
meansnorm = means / means[:,0].reshape((-1,1))
#meansnorm[meansnorm>1] = np.nan
areasnorm = areas / areas[:,0].reshape((-1,1))
#areasnorm[areasnorm>1] = np.nan
vols = areas * means
volsnorm = (vols - vols[:,0].reshape((-1,1)) ) / vols[:,0].reshape((-1,1))
In [111]:
for i, item in enumerate(areasnorm):
plt.plot(np.arange(20,0,-1), item, "b-", alpha=0.005)
In [152]:
matplotlib.rcParams.update({'font.size': 7})
plt.figure(figsize=(14,10))
for i, item in enumerate(range(1300,1400)):
ax1 = plt.subplot(10,10,i+1)
plt.plot(np.arange(20,0,-1), means[item], "b-")
plt.ylim(0,30)
plt.grid()
plt.title(df.dtime.iloc[item] + ", " + str(df.label.iloc[item]), fontsize=7)
ax2 = ax1.twinx()
plt.semilogy(np.arange(20,0,-1), areas[item], "r-")
plt.ylim(10,10000)
plt.tight_layout()
In [174]:
from scipy.signal import argrelextrema
In [176]:
np.gradient(areas[item])
Out[176]:
In [177]:
argrelextrema(np.gradient(areas[item]), np.greater)[0]
Out[177]:
In [219]:
blacklist_hours = ["2016-06-29 02:00:00"
"2016-06-29 13:00:00",
"2016-06-29 14:00:00",
"2016-07-05 05:00:00",
"2016-07-05 16:00:00",
"2016-07-05 17:00:00"]
blacklist_days = ["2016-06-16", "2016-06-29", "2016-07-04", "2016-07-05"]
In [221]:
for day in blacklist_days:
df.mean_intensity.loc[day] = -9999
bigx = np.argsort(df.mean_intensity)[::-1]
plt.figure()
plt.plot( np.array(df.mean_intensity)[bigx])
plt.ylim(0,60)
Out[221]:
In [225]:
matplotlib.rcParams.update({'font.size': 7})
plt.figure(figsize=(14,10))
# Look at the 100 most intensive objects
for i, item in enumerate(bigx[0:100]):
ax1 = plt.subplot(10,10,i+1)
plt.plot(np.arange(20,0,-1), np.gradient(areas[item]), "b-")
# for local maxima
extr = argrelextrema(np.gradient(areas[item]), np.greater)[0]
plt.plot(np.arange(20,0,-1)[extr], np.gradient(areas[item])[extr], "bo")
plt.grid()
plt.ylim(0,20)
plt.title(df.dtime.iloc[item] + ", " + str(df.label.iloc[item]), fontsize=7)
ax2 = ax1.twinx()
plt.plot(np.arange(20,0,-1), areasnorm[item], "g-")
plt.grid()
#plt.plot(np.arange(20,0,-1), areas[item], "r-")
plt.ylim(1,10)
plt.tight_layout()
In [148]:
plt.figure(figsize=(6,6))
ax1 = plt.subplot(2,1,1)
ax2 = plt.subplot(2,1,2)
for item in range(len(means[1:])):
ax1.plot(np.arange(20,0,-1), meansnorm[item], "b-", alpha=0.005)
ax1.grid()
ax2.semilogy(np.arange(20,0,-1), areasnorm[item], "r-", alpha=0.005)
In [124]:
plt.figure(figsize=(6,6))
for i, item in enumerate(volsnorm):
plt.plot(np.arange(20,0,-1), item, "b-", alpha=0.005)
In [228]:
time_window = np.arange(1, 7)
tdeltas = [dt.timedelta(seconds=i*3600.) for i in time_window]
In [234]:
tmeans = np.zeros( (len(df), len(time_window)+1) )
tareas = np.zeros( (len(df), len(time_window)+1) )
_dtime = None
with h5py.File(hourlyfile, 'r') as f:
for i in range(len(df)):
dtime = dt.datetime.strptime(df.dtime.iloc[i], "%Y-%m-%d %H:%M:%S")
bottom, left, top, right = df.bbox.iloc[i][0], df.bbox.iloc[i][1], df.bbox.iloc[i][2], df.bbox.iloc[i][3]
tmeans[i, 0] = df.mean_intensity.iloc[i]
tareas[i, 0] = df.area.iloc[i]
if dtime != _dtime:
print("")
print(dtime, end="")
# Only process new hourly set for a new datetime
thewindow = [dtime + item for item in tdeltas]
daystrings = [item.strftime("%Y/%m/%d") for item in thewindow]
hours = [item.hour for item in thewindow]
hsets = np.zeros((len(time_window), 900, 900)) * np.nan
for i in range(len(time_window)):
try:
hsets[i] = f[dtime.strftime("%Y/%m/%d")][dtime.hour]
except KeyError:
continue
hsets = np.cumsum(hsets, axis=0)
threshregions = [get_regions(hset, 20.) for hset in hsets]
else:
print(".", end="")
for trix, tr in enumerate(threshregions):
for r in tr:
# Looking for region that contains core region
if (left >= r.bbox[1]) and \
(right <= r.bbox[3]) and \
(bottom >= r.bbox[0]) and \
(top <= r.bbox[2]):
# Found
tmeans[i,trix+1] = r.mean_intensity
tareas[i,trix+1] = r.area
_dtime = dtime
In [236]:
np.save("hdf/tmeans_2016", tmeans)
np.save("hdf/tareas_2016", tareas)
In [237]:
tmeansnorm = tmeans / tmeans[:,0].reshape((-1,1))
#meansnorm[meansnorm>1] = np.nan
tareasnorm = tareas / tareas[:,0].reshape((-1,1))
In [239]:
matplotlib.rcParams.update({'font.size': 7})
plt.figure(figsize=(14,10))
for i, item in enumerate(range(1300,1400)):
ax1 = plt.subplot(10,10,i+1)
plt.plot(np.arange(0,7), tmeans[item], "b-")
plt.grid()
plt.title(df.dtime.iloc[item] + ", " + str(df.label.iloc[item]), fontsize=7)
ax2 = ax1.twinx()
plt.plot(np.arange(0,7), tareas[item], "r-")
plt.tight_layout()
In [60]:
radolan_grid_xy = wradlib.georef.get_radolan_grid(900,900)
x = radolan_grid_xy[:,:,0]
y = radolan_grid_xy[:,:,1]
In [242]:
dtime = "2016-06-07 21:00:00"
dtime_ = dt.datetime.strptime(dtime, "%Y-%m-%d %H:%M:%S")
with h5py.File(hourlyfile, 'r') as f:
hset = f[dtime_.strftime("%Y/%m/%d")][dtime_.hour]
sub = df.loc[dtime]
cmap=plt.cm.nipy_spectral
norm = matplotlib.colors.BoundaryNorm(np.arange(0,21), cmap.N)
plt.figure(figsize=(8,8))
pm = plt.pcolormesh(np.ma.masked_array(hset, ~np.isfinite(hset)), cmap=cmap, norm=norm)
plt.xlabel("RADOLAN easting (km)")
plt.ylabel("RADOLAN northing (km)")
plt.colorbar(pm)
ax = plt.gca()
patches = []
for i in range(0,len(sub)):
polygon = Rectangle(
(sub.iloc[i]["bbox"][1], sub.iloc[i]["bbox"][0]), # (x,y)
sub.iloc[i]["bbox"][3]-sub.iloc[i]["bbox"][1], # height
sub.iloc[i]["bbox"][2]-sub.iloc[i]["bbox"][0] # width
)
patches.append(polygon)
p = PatchCollection(patches, facecolor="None", edgecolor="white", linewidth=2)
ax.add_collection(p)
for i in range(0,len(sub)):
plt.text(sub.iloc[i].centroid[1], sub.iloc[i].centroid[0], str(sub.iloc[i].label), color="red", fontsize=18)
In [132]:
sub.centroid
Out[132]:
In [38]:
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111, aspect="equal")
patches = []
for i in range(1,len(df)):
polygon = Rectangle(
df.iloc[i]["bbox"][0:2], # (x,y)
df.iloc[i]["bbox"][2]-df.iloc[i]["bbox"][0], # width
df.iloc[i]["bbox"][3]-df.iloc[i]["bbox"][1], # height
)
patches.append(polygon)
colors = 100*np.random.rand(len(patches))
p = PatchCollection(patches, alpha=0.4)
p.set_array(np.array(colors))
ax.add_collection(p)
plt.xlim(0,900)
plt.ylim(0,900)
#plt.draw()
Out[38]:
In [ ]:
proj = wradlib.georef.create_osr("dwd-radolan")
watersheds_shp = r"E:\src\git\heisterm_bitbucket\tsms_data\tsms-data-misc\shapefiles\watersheds_kocher.shp"
dataset, inLayer = wradlib.io.open_vector(watersheds_shp)
cats, ids = wradlib.georef.get_vector_coordinates(inLayer, dest_srs=proj,
key="value")
ids = np.array(ids)
left, right, bottom, top = inLayer.GetExtent()
In [96]:
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111, aspect="equal")
patches = []
for i in range(1,len(df)):
polygon = Rectangle(
df.iloc[i]["bbox"][0:2], # (x,y)
df.iloc[i]["bbox"][2]-df.iloc[i]["bbox"][0], # width
df.iloc[i]["bbox"][3]-df.iloc[i]["bbox"][1], # height
)
patches.append(polygon)
colors = 100*np.random.rand(len(patches))
p = PatchCollection(patches, alpha=0.4)
p.set_array(np.array(colors))
ax.add_collection(p)
#for i in range(1,len(df)):
# plt.plot(df.ix[i]["centroid"][0], df.ix[i]["centroid"][1], "bo")
plt.xlim(0,900)
plt.ylim(0,900)
#plt.draw()
#wradlib.vis.add_lines(ax, cats, color='red', lw=0.5, zorder=4, alpha=0.3)
#plt.xlim(-40,20)
#plt.ylim(-4440,-4390)
Out[96]:
In [84]:
toobigx = np.where(df["area"]>1500)[0]
print(len(toobigx))
In [85]:
for i in toobigx:
plt.figure()
plt.pcolormesh(df.iloc[i]["image"])
In [77]:
i
Out[77]:
In [155]:
plt.hist(df["area"], bins=100, range=(0,200), log=True)
Out[155]:
In [ ]:
plt.figure(figsize=(8,8))
#plt.imshow(im.mean(axis=0), cmap=plt.cm.gray, origin="lower")
plt.imshow(frames[start:end].sum(axis=0), cmap=plt.cm.gray, origin="lower", vmin=0, vmax=30)
plt.xlabel("RADOLAN easting (km)")
plt.ylabel("RADOLAN northing (km)")
plt.title("Rainfall accumulation and cell tracks\nMay 29, 2016, 15:00-18:00 UTC")
ax = plt.gca()
for label in labels[1:]:
#for i in range(len(im)):
tmp = (label_im == label).astype("int")
#tmp = label_im[i]
regions = regionprops(tmp, intensity_image=im)
centx, centy = [], []
for region in regions:
y0, x0 = region.centroid
centx.append(x0)
centy.append(y0)
orientation = region.orientation
angle=-np.rad2deg( orientation)
e = Ellipse([x0,y0], region.major_axis_length, region.minor_axis_length,
angle=angle, facecolor="none", edgecolor="blue", linewidth=1.3, alpha=0.5)
ax.add_artist(e)
#plt.plot(x0, y0, "o", markerfacecolor=plt.cm.rainbow(i/len(im)), markeredgecolor="none", alpha=0.5)
plt.contour(tmp, [0.5], linewidths=1., colors="red", alpha=0.5)
#pm=plt.scatter([], [], c=[], cmap=plt.cm.rainbow, vmin=0, vmax=len(im)*5)
#cb=plt.colorbar(pm, label="Minutes from 2016-05-29 16:00", shrink=0.75)
In [ ]:
fig = plt.figure(figsize=(12,12))
ax = fig.add_subplot(111)
for i, label in enumerate(labels):
tmp = (label_im == label)
areal_avg = np.array([np.mean(frames[i][tmp]) for i in range(len(frames))])
ax.plot(np.cumsum(areal_avg))
In [ ]:
from matplotlib import animation
In [ ]:
# Animate features
# Prepare canvas
fig = plt.figure(figsize=(5,5))
ax = plt.subplot(111,aspect="equal")
im1 = ax.imshow(frames[0], origin="lower", cmap="gray", interpolation="none", vmin=10, vmax=20)
plt.xlabel("Easting (km)")
plt.ylabel("Northing (km)")
plt.grid(color="white")
plt.xlim(150,450)
plt.ylim(550,900)
#ax1.plot(x[0,goodtrack], y[0,goodtrack], linestyle="None", marker="o", mfc="None", mec="limegreen")
#ax1.plot(x[0,~goodtrack], y[0,~goodtrack], linestyle="None", marker="o", mfc="None", mec="red")
ax.grid(color="white")
tstamp1 = ax.text(160, 560, dtimes[0].isoformat(), color="white", fontsize=12)
def animate(j):
im1.set_array(frames[0+j])
tstamp1.set_text(dtimes[0+j].isoformat())
for label in labels[1:]:
#break
tmp = (label_im[j] == label).astype("int")
#tmp = label_im[i]
regions = regionprops(tmp, intensity_image=im[j])
centx, centy = [], []
for region in regions:
y0, x0 = region.centroid
centx.append(x0)
centy.append(y0)
orientation = region.orientation
angle=-np.rad2deg( orientation)
e = Ellipse([x0,y0], region.major_axis_length, region.minor_axis_length,
angle=angle, facecolor="none", edgecolor=plt.cm.rainbow(j/len(im)), linewidth=1.3, alpha=0.3)
ax.add_artist(e)
#ax.plot(x0, y0, "o", markerfacecolor=plt.cm.rainbow(j/len(im)), markeredgecolor="none", alpha=0.5)
tstamp1.set_text(dtimes[0+j].isoformat())
return im1
# ATTENTION: THIS IS SLOW - Rendering each frame of the animation might take more time than the interval between the frames
# This can cause the temporal sequence to be confused in the matplotlib interactive mode.
# The animation thus looks better if saved as movie, or you have to increase the interval argument
# Animation not shown in notebook if you use %pylab inline
maxi = len(frames)-1
ani = animation.FuncAnimation(fig, animate, frames=np.arange(0, maxi), interval=400, blit=False)
ani.save("features.gif", writer="imagemagick", dpi=150)
In [ ]:
len(region)
In [ ]:
#fig, ax = plt.subplots()
plt.imshow(im, cmap=plt.cm.gray, origin="lower")
plt.contour(label_im, [0.5], linewidths=1.2, colors='y')
plt.xlabel("RADOLAN easting (km)")
plt.ylabel("RADOLAN northing (km)")
plt.title("Snaphot at 2016-05-29 16:00 UTC")
ax = plt.gca()
for i, props in enumerate(regions):
y0, x0 = props.centroid
orientation = props.orientation
x1 = x0 + math.cos(orientation) * 0.5 * props.major_axis_length
y1 = y0 - math.sin(orientation) * 0.5 * props.major_axis_length
x2 = x0 - math.sin(orientation) * 0.5 * props.minor_axis_length
y2 = y0 - math.cos(orientation) * 0.5 * props.minor_axis_length
#plt.plot((x0, x1), (y0, y1), '--r', linewidth=2)
#plt.plot((x0, x2), (y0, y2), '--r', linewidth=2)
#plt.plot(x0, y0, '.r', markersize=15)
angle=-np.rad2deg( props.orientation)
e = Ellipse([x0,y0], props.major_axis_length, props.minor_axis_length,
angle=angle, facecolor="none", edgecolor="red", linewidth=2)
ax.add_artist(e)
minr, minc, maxr, maxc = props.bbox
bx = (minc, maxc, maxc, minc, minc)
by = (minr, minr, maxr, maxr, minr)
#plt.plot(bx, by, '-b', linewidth=2.5)
try:
label = "ID=%s\navg=%d mm/h\nmax=%d mm/h" % (props.label, props.mean_intensity, props.max_intensity)
except:
label = "ID=%s, avg=%s mm/h, max=%s mm/h" % (props.label, "nan", "nan")
plt.text((minc+maxc)/2, maxr+2, label, color="red", fontsize=10, horizontalalignment='center')
#plt.axis((0, 900, 900, 0))
plt.xlim(200,900)
plt.ylim(0,470)
In [ ]:
minr, minc, maxr, maxc = props.bbox
plt.imshow(im[minr:maxr, minc:maxc])
In [ ]:
im2 = frames[1]
fig = plt.figure(figsize=(8, 8))
ax2 = plt.subplot(1, 1, 1)
for i, props in enumerate(regions):
minr, minc, maxr, maxc = props.bbox
roi = im[minr:maxr, minc:maxc]
result = match_template(im2, roi)
ij = np.unravel_index(np.argmax(result), result.shape)
x, y = ij[::-1]
print(ij)
#ax1.imshow(roi, cmap=plt.cm.gray)
#ax1.set_axis_off()
#ax1.set_title('Feature #1 at t+0')
ax2.imshow(im2, cmap=plt.cm.gray, origin="lower")
ax2.set_axis_off()
ax2.set_title('Feature #1 at t+2')
# highlight matched region
hcoin, wcoin = roi.shape
rect = plt.Rectangle((x, y), wcoin, hcoin, edgecolor='r', facecolor='none')
ax2.add_patch(rect)
plt.plot(x,y,".r")
plt.plot(ij[0],ij[1],".b")
# highlight matched region
bx = (minc, maxc, maxc, minc, minc)
by = (minr, minr, maxr, maxr, minr)
plt.plot(bx, by, '-b', linewidth=1.)
In [ ]:
ij
In [ ]:
ndimage.find_objects(label_im==15)
In [ ]:
image = frames[2]
coin = roi
result = match_template(image, coin)
ij = np.unravel_index(np.argmax(result), result.shape)
x, y = ij[::-1]
fig = plt.figure(figsize=(8, 3))
ax1 = plt.subplot(1, 2, 1)
ax2 = plt.subplot(1, 2, 2, adjustable='box-forced')
ax1.imshow(coin, cmap=plt.cm.gray)
ax1.set_axis_off()
ax1.set_title('Feature #1 at t+0')
ax2.imshow(image, cmap=plt.cm.gray)
ax2.set_axis_off()
ax2.set_title('Feature #1 at t+2')
# highlight matched region
hcoin, wcoin = coin.shape
rect = plt.Rectangle((x, y), wcoin, hcoin, edgecolor='r', facecolor='none')
ax2.add_patch(rect)
In [ ]: