In [1]:
%matplotlib
import wradlib
#import cv2
import numpy as np
import os
import matplotlib.pyplot as plt
#from matplotlib import animation
import matplotlib.patches as mpatches
from matplotlib.collections import PatchCollection
from scipy.ndimage import zoom
import datetime
import warnings
warnings.simplefilter('once', DeprecationWarning)
In [2]:
from scipy import ndimage as ndi
from skimage import feature
from skimage.feature import match_template
In [3]:
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
Data is from the German Weather Service: the so called RY product represents rainfall intensity composite for the whole of Germany in 5 minute intervals.
Spatial resolution: 1 x 1 km; spatial extent: 900 x 900 km.
Information required from user
datadir where you store the RY data (unpack the ry archives there).dtimes lines.downsizeby) in order to avoid memory problems (this becomes relevant once you solve the 2D-adveciton equation...)
In [29]:
# Set data directory
datadir = "data/ry"
# Original grid dimensions
nx = 900
ny = 900
# pixel size (in meters)
dx = 1000.
dy = 1000.
# Downsize by factor "downsizeby"
# downsizeby = 1 will leave the dimensions unchanged,
# but for a 900x900 km grid, downsizing might be
# required in order to avoid MemoryError
downsizeby = 1
# interval between observations (in seconds)
interval = 300
# Set time window
##dtimes = wradlib.util.from_to("2008-06-02 17:00:00", "2008-06-02 19:00:00", interval)
##dtimes = wradlib.util.from_to("2015-04-26 17:00:00", "2015-04-26 19:00:00", interval)
##dtimes = wradlib.util.from_to("2015-03-29 17:00:00", "2015-03-29 19:00:00", interval)
##dtimes = wradlib.util.from_to("2016-05-29 16:30:00", "2016-05-29 18:30:00", interval)
dtimes = wradlib.util.from_to("2016-05-23 04:00:00", "2016-05-23 08:00:00", interval)
In [30]:
# Compute grid dimensions and grid coordinates after resampling
dx2, dy2 = dx*downsizeby, dy*downsizeby
nx2, ny2 = int(nx/downsizeby), int(ny/downsizeby)
X2, Y2 = np.meshgrid( np.arange(0,nx2*dx2, dx2), np.arange(0,ny2*dy2, dy2) )
# Define container
frames = np.zeros( (len(dtimes), nx2, ny2 ) )
# Read the data, convert back to dBZ, and downsize
# (maybe also try with keeping mm/h instead of converting to dBZ?)
for i, dtime in enumerate(dtimes):
fname = dtime.strftime( os.path.join(datadir, "raa01-ry_10000-%y%m%d%H%M-dwd---bin") )
frames[i] = zoom( wradlib.io.read_RADOLAN_composite(fname, missing=0)[0], 1./downsizeby, order=1)
frames[i] = wradlib.trafo.decibel( wradlib.zr.r2z(frames[i]) )
frames[i][frames[i]<0] = 0
In [31]:
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_shape(watersheds_shp)
cats, ids = wradlib.georef.get_shape_coordinates(inLayer, dest_srs=proj,
key="value")
ids = np.array(ids)
left, right, bottom, top = inLayer.GetExtent()
In [32]:
radolan_grid_xy = wradlib.georef.get_radolan_grid(900,900)
x = radolan_grid_xy[:,:,0]
y = radolan_grid_xy[:,:,1]
In [1]:
fig = plt.figure(figsize=(12, 12))
ax = fig.add_subplot(111)
plt.pcolormesh(x, y, frames.sum(axis=0), vmin=0, cmap="gray")
cb = plt.colorbar(shrink=0.75)
wradlib.vis.add_lines(ax, cats, color='red', lw=0.5, zorder=4, alpha=0.3)
#plt.xlim(-40,20)
#plt.ylim(-4440,-4390)
In [33]:
#im = wradlib.trafo.r2depth(frames[12:12+24].mean(axis=0), interval=300*20)
im = frames
In [34]:
blobs = im > 15.
label_im = measure.label(blobs, background=0)
nb_labels = len(np.unique(label_im))
#sizes = ndimage.sum(blobs, label_im, range(nb_labels + 1))
#mask_size = sizes < 1500
#remove_pixel = mask_size[label_im]
#label_im[remove_pixel] = 0
#labels = np.unique(label_im)
#label_im = np.searchsorted(labels, label_im)
regions = regionprops(label_im, intensity_image=im)
print(len(regions))
In [35]:
minsize = 100
maxsize = 4000
minintensity = 22
labels = np.array([region.label for region in regions])
islarge = np.array([(region.area >= minsize) for region in regions])
#islarge = np.array([((region.area >= minsize) & (region.area < maxsize)) for region in regions])
largelabels = labels[islarge]
print(len(largelabels))
In [36]:
label_im[~np.isin(label_im, largelabels)] = 0
labels = np.unique(label_im)
In [12]:
plt.figure(figsize=(8,8))
#plt.imshow(im.mean(axis=0), cmap=plt.cm.gray, origin="lower")
plt.imshow(frames.mean(axis=0), cmap=plt.cm.gray, origin="lower", vmin=15)
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[i] == label).astype("int")
#tmp = label_im[i]
regions = regionprops(tmp, intensity_image=im[i])
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(i/len(im)), 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=[plt.cm.spectral(i/len(im)),], 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 [13]:
from matplotlib import animation
In [38]:
# 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 [ ]: