(1) If the eigenvalues of $f$ always have multiplicity 1, then indeed you can make a continuous selection of eigenvalueseigenvectors. Say $\lambda_j(z)$, with $j=1,2,\dots,n$, denote the eigenvalues of $f(z)$ arranged in increasing order. For each $j$ choose a continuous $\alpha_j(t,z)$, with $t\in \mathbb R$ and $z\in D$ such that $$\alpha_j(\lambda_j(z),z)=1\mbox{ and }\alpha_j(\lambda_{i}(z),z)=0\mbox{ if }i\neq j.$$ Using functional calculus let us define $$p_j(z)=\alpha_j(f_j(z),z).$$ These are the rank one orthogonal projections onto the eigenspaces of $f(z)$. Since $D$ is an open subset of the plane, it only has trivial complex line bundles. So for each each $p_j(z)$ there exists a continuous section $v_j(z)\in \mathbb C^n$ such that $p_j(z)v_j(z)=v_j(z)$.
(2) The set of functions $f$ with all eigenvalues of multiplicity 1 is dense (and $G_\delta$) among the bounded continuous functions on $D$ with values on the $n\times n$ selfadjoint matrices. Here it is crucial that $D$ is at most of dimension 2. This is proven in "Density of the self-adjoint elements with finite spectrum in an irrational rotation Cā-algebra. Math. Scand. 67 (1990), 73ā86.", by Choi and Elliott.
The gist of their argument is this: the set of $n\times n$ self-adjoint matrices such that at least two eigenvalues agree is a finite union of submanifolds of the set of all self-adjoint matrices, where each submanifold has codimension at least 3. This makes it that a suitable perturbation of $f$ can avoid such finite union of submanifolds provided that the domain of $f$ has dimension at most 2.
