The Riemann integral of $f$ over an arbitrary bounded region $D \subset \mathbb{R}^2$ is defined as
$$\int_Df = \int_Q\hat{f},$$
$\hat{f}(x) = f(x)$ for $x \in D$ and $\hat{f}(x) = 0$ for $x \notin D$, and $Q$ is any rectangle containing $D$. It is straightforward to show that the choice for $Q$ is arbitrary.
Assume that $\hat{f}$ is Riemann integrable over $Q$. This means there exists a real number $I$ -- the integral-- such that for any $\epsilon > 0$ there exists a partition $P$ with subrectangles $\{R_1, \ldots, R_n\}$ such that
$$\tag{1}- \frac{\epsilon}{2} < S(P,f) - I < \frac{\epsilon}{2}$$
Here $S(P,f) = \sum_{j=1}^n f(t_j) \, vol(R_j)$ is a Riemann sum with any choice of intermediate points $t_j \in R_j$ for $j = 1, \ldots, n$.
Recall that upper and lower Darboux sums are given by
$$U(P,f) = \sum_{j=1}^n M_j\, vol(R_j), \quad M_j := \sup_{x \in R_j}f(x),\\ L(P,f) = \sum_{j=1}^n m_j\, vol(R_j), \quad m_j := \inf_{x \in R_j}f(x)$$
By the properties of supremum and infimum there exists $\xi_j, \eta_j \in R_j$ such that
$$M_j - \frac{\epsilon}{2 \,vol(Q)} < f(\xi_j) \leqslant M_j, \quad m_j \leqslant f(\xi_j) < m_j + \frac{\epsilon}{2 \,vol(Q)}$$
Summing over all subrectangles of the partition we get
$$\tag{2}U(P,f) - \frac{\epsilon}{2} < \sum_{j=1}f(\xi_j) \, vol(R_j), \quad \sum_{j=1}f(\eta_j) \, vol(R_j)< L(P,f) + \frac{\epsilon}{2}$$
The sums in (2) are Riemann sums with intermediate points $\xi_j$ and $\eta_j$, respectively. Applying the inequalities in (1), it follows that $U(P,f) - \frac{\epsilon}{2} < I + \frac{\epsilon}{2}$ and $I - \frac{\epsilon}{2} < L(P,f) - \frac{\epsilon}{2}$.
Hence,
$$I - \epsilon < L(P,f) \leqslant U(P,f) < I + \epsilon$$
There are several ways to proceed in showing that this proves Darboux integrability. One is to note that the upper and lower Darboux integrals are sandwiched between any upper and lower Darboux sums, and, thus,
$$I- \epsilon < \underline{\int_Q} \hat{f} \leqslant \overline{\int_Q} \hat{f} < I + \epsilon$$
Since $\epsilon > 0$ can be chosen arbitrarily close to $0$, Darboux integrability holds with
$$\underline{\int_Q} \hat{f} = \overline{\int_Q} \hat{f}$$
Best Answer
Let $\mathbf E$ be a Banach space, let $a<b$ be real numbers, let $f : [a,b] \to \mathbf E$ be a function.
A partition $\pi$ of $[a,b]$ is a finite subset, $\{a,b\} \subseteq \pi \subset [a,b]$, usually written in order: $a = x_0 < x_1 < \dots < x_n = b$. Tags for a partition $\pi$ as above are points $t_i$ such that $x_{i-1} \le t_i \le x_i$ for $1 \le i \le n$. Partition $\pi_1$ refines partition $\pi_2$ iff $\pi_1 \supseteq \pi_2$, remembering that a partition is a finite set.
Definition. Let $f$ be as above, and let $\mathbf u \in \mathbf E$. We say that $f$ is Riemann integrable on $[a,b]$ and $\mathbf{u}$ is its integral iff: for every $\epsilon > 0$, there is a partition $\pi_0$ of $[a,b]$ such that for all refinements $\pi=(x_i)_{i=0}^n$ of $\pi_0$ and all tags $(t_i)_{i=1}^n$ for $\pi$, $$ \left\|\mathbf u - \sum_{i=1}^n f(t_i)\;(x_{i}-x_{i-1})\right\| < \epsilon. $$
Lemma $f$ is integrable iff: for every $\epsilon > 0$ there is a partition $\pi = (x_i)_{i=0}^n$ such that for any two choices $(t_i)_{i=1}^n, (s_i)_{i=1}^n$ of tags for $\pi$, we have $$ \left\|\sum_{i=1}^n \big(f(t_i)-f(s_i)\big)\;(x_{i}-x_{i-1})\right\| < \epsilon. $$
Proof. Cauchy criterion.
Theorem. Let $f : [a,b] \to \mathbf E$ be bounded and continuous except on a set $N\subseteq [a,b]$ of measure zero. Then $f$ is integrable.
Proof. Add $\{a,b\}$ to the null set $N$ to avoid special cases for endpoints. Let $\epsilon>0$. Say $f$ is bounded by $M$, $\|f(x)\| \le M$. Let $\alpha > 0$ be so small that $2M\alpha + \alpha(b-a) < \epsilon$. For an open interval $(u,v)$ we say $f$ has oscillation at most $\alpha$ on $(u,v)$ if for all $x,y \in (u,v)$, $\|f(x)-f(y)\| \le \alpha$. If $f$ is continuous at a point $s$, then there is an inverval $(u,v)$ with rational endpoints, $s \in (u,v)$, so that $f$ has oscillation at most $\alpha$ on $(u,v)$. So there is a countable union of such intervals $(u,v)$ that contains $[a,b] \setminus N$, and thus has full measure. So there is a finite list $(u_j,v_j)$ of intervals where $f$ has oscillation at most $\alpha$, and their union has measure greater than $b-a-\alpha$. Then there is a partition $\pi = (x_i)_{i=0}^n$ of $[a,b]$ such that each subinterval $[x_{i-1},x_i]$ from the partition either is contained in an interval where $f$ has oscillation at most $\alpha$, or is an "exceptional" interval. The total length of all the exceptional intervals is ${}< \alpha$. Now let $(t_i)$ and $(s_i)$ be two choices of tags for the partition $\pi$. Now we must consider $$ \left\|\sum_{i=1}^n\big(f(t_i)-f(s_i)\big)(x_i-x_{i-1})\right\| \le \sum_{i=1}^n\left\|\big(f(t_i)-f(s_i)\big)(x_i-x_{i-1})\right\| . $$ Consider the term $(f(t_i)-f(s_i))(x_i-x_{i-1})$. If the subinterval $[x_{i-1},x_i]$ is not exceptional, then $$ \left\|\big(f(t_i)-f(s_i)\big)(x_i-x_{i-1})\right\| \le \alpha (x_{i}-x_{i-1}) , $$ so the total of all terms for non-exceptional intervals is at most $\alpha (b-a)$. If the subinterval $[x_{i-1},x_i]$ is exceptional, then $$ \left\|\big(f(t_i)-f(s_i)\big)(x_i-x_{i-1})\right\| \le 2M (x_{i}-x_{i-1}) , $$ so the total of all terms for exceptional intervals is at most $2M\alpha$. Thus $$ \left\|\sum_{i=1}^n\big(f(t_i)-f(s_i)\big)(x_i-x_{i-1})\right\| \le \alpha(b-a)+2M\alpha < \epsilon . $$