I just happen to have written code to do exactly this yesterday.
Think in terms of subtracting sets, and think recursively.
For sets $A_1 ... A_n$, compute $A_1' = A_1$, then $A_2' = A_2 \backslash A_1'$, then $A_3' = (A_3 \backslash A_1') \backslash A_2'$, and so on. By a simple inductive argument, $A_1'...A_n'$ are disjoint. If it helps, here's a Haskell function that does this:
disjoint [] = []
disjoint (a:as) = a : disjoint (concatMap (\\ a) as)
where "\\" is a set subtraction operator that returns a list of disjoint sets (to represent a union), and (by Haskell's sectioning rules) "(\\ a)" is a function that subtracts "a" from its input. This function preserves the head of the list of sets, subtracts the head from the rest of them, flattens the resulting list of lists, and applies itself to the result.
For $d$-dimensional sets, think recursively again. A $d$-dimensional rectangle is a cartesian product $I_d \times \prod_{j=1}^{d-1} I_j$, where $I_d$ is, say, a one-dimensional thing like an interval. Thus, you need functions for subtracting intervals and for subtracting two-dimensional products of any kind of set, such as products of two intervals, products of an interval and a product, etc.
Subtracting intervals isn't too hard, but it's a bit tedious to work out all the cases if you track whether endpoints are open or closed.
You can turn the following formula into a function for subtracting products:
$$(I_1 \times I_2) \backslash (J_1 \times J_2) = ((I_1 \backslash J_1) \times J_2) \uplus ((I_1 \cap J_1) \times (I_2 \backslash J_2))$$
Because subtraction can yield a list (i.e. union) of intervals or products, it's a little trickier to implement than just typing it in, but it's not too hard. Haskell code again, in case there are native speakers about:
PairSet a1 a2 \\ PairSet b1 b2 =
let cs = do c1 <- a1 \\ b1
return (PairSet c1 a2)
ds = do d2 <- a2 \\ b2
return (PairSet (a1 /\ b1) d2)
in cs ++ ds
where "/\" computes intersections of arbitrary sets.
So if you have subtraction and intersection (which can be defined in terms of subtraction) for intervals and products, something like the above code will recursively subtract rectangular sets of any data shape you have, with $d$-dimensional rectangles (i.e. lists of intervals) as a special case.
Now, "disjoint" is obviously $O(n^2)$, which isn't bad. But each interval subtraction can generate two intervals, so you're looking at exponential time (in $d$) in the worst case. I haven't discovered any worst cases yet, but I haven't tried very hard yet.
The code above distributes lists (i.e. unions) of intervals over products. In some past code, I kept this from happening by representing sets of reals as sorted, disjoint intervals instead of just as intervals, but I've never measured how much it helps. If it does, you could probably do even better by representing sets of pairs as a union of products in a BSP tree, as Joseph suggests. It would basically be a method to compress disjoint unions of rectangles, in a way that doesn't complicate common operations (i.e. product, union, intersect, subtract, measure volume).
