In [47]:

    

url = 'http://opendap-matroos.deltares.nl/thredds/catalog/maps2d/knmi_h11_v72/catalog.xml'
import osgeo.osr


    

In [48]:

    

from thredds_crawler.crawl import Crawl
import numpy as np
import matplotlib.pyplot as plt
import cmocean.cm
import matplotlib.colors
# %matplotlib inline


    

In [3]:

    

c = Crawl(url)


    

In [4]:

    

dataset = c.datasets[0]


    

In [5]:

    

service = [service for service in dataset.services if service['service'] == 'OPENDAP'][0]


    

In [49]:

    

url = service['url']
url
url = '/Users/baart_f/data/knmi/201509230600.nc'


    

In [50]:

    

import netCDF4


    

In [51]:

    

ds = netCDF4.Dataset(url)


    

In [52]:

    

u = ds.variables['northward_wind']
v = ds.variables['eastward_wind']


    

In [ ]:

    




    

In [53]:

    

def plot(t):
    uv = np.sqrt(u[t]**2 + v[t]**2)
    plt.imshow(uv, cmap=cmocean.cm.option_d)
    plt.colorbar()
from IPython.html.widgets import interactive
interactive(plot, t=(0, u.shape[0]))


    

In [54]:

    

X = ds.variables['x'][:,:]
Y = ds.variables['y'][:,:]


    

In [55]:

    

# for i in range(u.shape[0]):
i = 0
uv = np.sqrt(u[i]**2 + v[i]**2)
plt.figure(figsize=(13,8))
# cmocean.cm.make_speed_cmap()
plt.pcolormesh(X, Y, uv, cmap=cmocean.cm.option_d, vmin=0, vmax=20, shading='gouraud')
plt.axis('off')
plt.tight_layout()
plt.savefig('knmi_%04d.png' % (i,))


    

In [56]:

    

# let's define the systems
src_srs = osgeo.osr.SpatialReference()
src_srs.ImportFromEPSG(4326)
dst_srs = osgeo.osr.SpatialReference()
dst_srs.ImportFromEPSG(900913)
wgs84 = osgeo.osr.SpatialReference()
wgs84.ImportFromEPSG(4326)
utm = osgeo.osr.SpatialReference()
utm.ImportFromEPSG(32631)

# and the translations between them
src2wgs84 = osgeo.osr.CoordinateTransformation(src_srs, wgs84)
dst2wgs84 = osgeo.osr.CoordinateTransformation(dst_srs, wgs84)
utm2wgs84 = osgeo.osr.CoordinateTransformation(utm, wgs84)
wgs842utm = osgeo.osr.CoordinateTransformation(wgs84, utm)
utm2dst = osgeo.osr.CoordinateTransformation(utm, dst_srs)
src2utm = osgeo.osr.CoordinateTransformation(src_srs, utm)

src2dst = osgeo.osr.CoordinateTransformation(src_srs, dst_srs)


    

In [57]:

    

points_gm = np.array(src2dst.TransformPoints(np.c_[X.ravel(), Y.ravel()]))


    

In [58]:

    

XYZ_gm = points_gm.reshape(X.shape + (3,))
X_gm, Y_gm = XYZ_gm[...,0], XYZ_gm[...,1]


    

In [59]:

    

U, V = u[-1], v[-1]
plt.hist(U.ravel())
plt.hist(V.ravel())


    

    
Out[59]:





(array([  1229.,   4530.,  11585.,  29333.,  62586.,  28273.,  14284.,
         10451.,   4154.,    896.]),
 array([-10.77807617,  -8.38201971,  -5.98596325,  -3.58990679,
         -1.19385033,   1.20220613,   3.5982626 ,   5.99431906,
          8.39037552,  10.78643198,  13.18248844]),
 <a list of 10 Patch objects>)

In [60]:

    

cmap = cmocean.cm.option_d
# normalization function


#plt.imshow(C)


    

In [61]:

    

#def uv2rgb(X, Y, uv):
def uv2c(U, V):
    N = matplotlib.colors.Normalize(vmin=-15, vmax=15, clip=True)
    C = np.dstack([
            N(U), 
            N(V), 
            np.zeros_like(V),
            np.ones_like(V)
        ]).astype('float64')
    #C = cmap(N(U))
    C.min(), C.max(), C.shape, C.dtype
    return C

def fig2rgba_array(fig):
    fig.canvas.draw()
    buf = fig.canvas.tostring_argb()
    ncols, nrows = fig.canvas.get_width_height()
    arr = np.fromstring(buf, dtype=np.uint8).reshape(nrows, ncols, 4)
    # argb -> rgba (roll back 1)
    arr = np.roll(arr, -1, -1)
    return arr

def uv2rgba(X, Y, U, V, filename, res=1024):
    C = uv2c(U, V)
    fig, ax = plt.subplots(figsize=(res/72.0, res/72.0))
    ax.axis('off')
    fig.tight_layout(pad=0)
    plt.subplots_adjust(left=0, right=1, top=1, bottom=0)

    im = ax.pcolormesh(X_gm, Y_gm, np.zeros_like(X),  shading='gouraud')

    im.set_array(None) # remove the array
    im.set_edgecolor('none') # no borders
    # Pass the image, make sure it is type: float, shape N,M,4 (alpha channel)
    im.set_facecolor(C.reshape(C.shape[0]*C.shape[1], C.shape[2]))

    xlim = ax.get_xlim()
    ylim = ax.get_ylim()
    plt.plot(xlim[0], ylim[0], 'k+', markersize=100)
    plt.plot(xlim[0], ylim[1], 'k+', markersize=100)
    plt.plot(xlim[1], ylim[1], 'k+', markersize=100)
    plt.plot(xlim[1], ylim[0], 'k+', markersize=100)
    ax.set_xlim(xlim)
    #flip axis for webgl
    ax.set_ylim(ylim[::-1])
    fig.savefig(filename,  facecolor='blue', dpi=72)
    plt.close(fig)

n = u.shape[0]
step = 5
for i in range((n - 1)*step -1):
    i0, remainder = divmod(i, step)
    i1 = i0 + 1
    # interpolation factors
    f0 = (1-remainder/float(step))
    f1 = (remainder/float(step))
    # interpolate in time
    U = u[i0] * f0 + u[i1] * f1
    V = v[i0] * f0 + v[i1] * f1    
    uv2rgba(X, Y, U, V, 'knmi_%04d.png' % (i,))


    

In [27]:

    

fig = plt.figure()
fig.savefig?


    

In [32]:

    

u.shape[0]


    

    
Out[32]:





49

In [37]:

    

ds.variables['time'][:]


    

    
Out[37]:





array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
       17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,
       34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48], dtype=int32)

In [46]:

    

[X.min(), Y.min()], [X.max(), Y.max()]


    

    
Out[46]:





([-30.868795, 26.064999], [64.929413, 80.0])

In [45]:

    

ds


    

    
Out[45]:





<type 'netCDF4.Dataset'>
root group (NETCDF3_CLASSIC data model, file format DAP2):
    history: Data retrieved from KNMI (NL); Created by Matroos on Wed Sep 23 10:44:53 CEST 2015
    grid_type: IRREGULAR
    coordinate_system: WGS84
    Conventions: CF-1.0
    DODS_EXTRA.Unlimited_Dimension: time
    dimensions(sizes): time(49), x(371), y(451)
    variables(dimensions): float32 sea_ice_area_fraction_mean_sealevel(y,x), int32 time(time), float32 surface_snow_amount(time,y,x), float32 surface_net_upward_shortwave_flux(time,y,x), float32 x(y,x), float32 air_pressure(time,y,x), float32 air_temperature_at_cloud_top(time,y,x), float32 surface_roughness_length_mean_sealevel(y,x), float32 precipitation_amount(time,y,x), float32 y(y,x), float32 atmosphere_boundary_layer_thickness(time,y,x), float32 dew_point_temperature(time,y,x), float32 land_area_fraction(y,x), float32 cloud_area_fraction(time,y,x), float32 surface_upward_sensible_heat_flux(time,y,x), float32 low_cloud_area_fraction(time,y,x), float32 medium_cloud_area_fraction(time,y,x), float32 g10_rot_3(y,x), float32 air_pressure_fixed_height(time,y,x), float32 surface_net_upward_longwave_flux(time,y,x), float32 surface_snow_thickness(time,y,x), float32 high_cloud_area_fraction(time,y,x), float32 geopotential(y,x), float32 surface_upward_latent_heat_flux(time,y,x), float32 convective_precipitation_amount(time,y,x), float32 surface_albedo(time,y,x), float32 surface_roughness_length(y,x), float32 surface_downwelling_shortwave_flux_in_air(time,y,x), float32 air_temperature(time,y,x), float32 northward_wind(time,y,x), float32 eastward_wind(time,y,x)
    groups: 

In [ ]: