Ndsolve of random (periodic) initial conditions via BSplineFunction

Suppose we want to solve the heat equation for a random initial temperature field T(t,x,y). Therefore I define the initial conditions with the BSplineFunction which nicely incorporates periodic boundary conditions and try to solve it with NDSolve:

Clear["Global`*"] Lx = 1; Ly = 1; eq = D[T[t, x, y], t] == Laplacian[T[t, x, y], {x, y}]; pbc1 = T[t, x, 0] == T[t, x, Ly]; pbc2 = T[t, 0, y] == T[t, Lx, y]; T0[x_, y_] =   BSplineFunction[RandomReal[1, {30, 30, 1}], SplineClosed -> True][x,    y] Plot3D[T0[x, y], {x, 0, Lx}, {y, 0, Ly}, PlotRange -> All,   AxesLabel -> {"x", "y", "T"}, ColorFunction -> "DarkRainbow",   Mesh -> All, MeshStyle -> Opacity[.2]] T0[0, 0] ic = T[0, x, y] == T0[x, y] ic2 = T[0, x, y] == 0.5*Exp[-4 ((x - Lx/2)^2 + (y - Ly/2)^2)]; Monitor[AbsoluteTiming[   sol = NDSolve[{eq, pbc1, pbc2, ic},      T, {t, 0, 1}, {x, 0, Lx}, {y, 0, Ly},      Method -> {"MethodOfLines",        "SpatialDiscretization" -> {"TensorProductGrid",          "MaxPoints" -> 25}},      EvaluationMonitor :> (currentTime =         Row[{"t = ", CForm[t]}])]], currentTime] 

The code unfortunately throws an error "Data … is not a rectangular tensor with dimensions …". I assume this results from the fact that T0 evaluated e.g. at grid point [0,0] returns a list {} and NDSolve probably needs just the float value.

I tried T0[0,0][[1]] which yields the expected value. But this does not work for the definition of the initial condition: ic = T[0, x, y] == T0[x, y][[1]], where x,y are not set yet. This returns just x, in reality it should first put in the values of x,y and only after that take list element 1 (if I am not mistaken).

Apart from that the code works fine which can be verified if you simply use ic2 instead of ic as initial condition in NDSolve.

Can someone help to get this running?