First of all, Kd-trees are data-structures for a number of points, while your problem is considered with shapes. Although you can (and we will) represent the shapes as points in a higher dimensional space in some way, querying shapes is an important special case over arbitrary points. I don't see why Z-orders would be interesting for points, so we really want a data-structure for shapes. R-trees store rectangles, but if you are interested in finding intersections of query regions with non-rectangular shapes, R-trees probably aren't very helpful.
Since a general shape is not so clearly enough, I will consider searching for a range in a set of $n$ simple polygons, this should be sufficient for most cases. Querying this structure with a 2D range (that is, an orthogonal rectangle) for shapes can mean two different things:
- Find all polygons fully contained in an orthogonal query rectangle.
- Find all polygons intersecting an orthogonal query rectangle.
Polygons fully contained in an orthogonal query rectangle
For the first query, note that a polygon is fully contained in a rectangle if and only if its rectangular bounding box is fully contained inside the query rectangle. So, we can search for the bounding boxes of the shapes. To do this, we represent the bounding box of every shape $s$ as a tuple of the bottom left and upper right corner $(p_s, q_s)$. I denote the $x$- and $y$-coordinates of a point using brackets, e.g. $p[x]$. We want to find all $s$ such that $p_r \leq p_s$ and $q_r \leq q_s$ for a query rectangle $r$. We can use a 4D version of your favorite point range query data structure (Kd,R,Range,etc.-tree) to find all those rectangles in $O(\log^3 n + k)$ time, were $k$ is the amount of rectangles reported, and use $O(n\log^3 n)$ storage for the data-structure. You can actually reduce the storage to $O(n\log^2 n)$, using a 2D-Range tree on the points $p_s$ and 'linking' this tree to a priority search tree for the $q_s$ coordinate (since we only have to search in positive direction on the second coordinate). However, Range trees tend to have large coefficients in practice, so you should probably stick to the trees you're familiar with.
Polygons intersecting an orthogonal query rectangle
For the second type of query on points, finding all polygons intersecting an orthogonal query rectangle, searching for intersecting bounding boxes may lead to a lot bounding boxes where the polygon does not actually intersect the query rectangle, so we need another representation. If we represent all polygons by the $m$ segments of their boundaries, we use a Windowing query to find the segments in $O(log^2 m + i)$, where $i$ is the amount of segments found, using a $O(m\log m)$ size data-structure.
However, we missed the polygons that contain our query rectangle, as we only looked for boundary intersections. We can check this by storing the bounding boxes of all polygons in 2 'chained' segment trees (similar to creating from 2D-range trees from BST's) and query for the bounding boxes that contain endpoints of the query rectangle. This has a query time of $O(\log n + k)$ and uses $O(n\log n )$ storage, but this is dominated by the other part, as $n\leq m$ and $k \leq i$.
Now, for z-orders...
First, to summarize, we now have algorithms to solve the following problems with the following bounds:
- Finding all polygons fully contained in an orthogonal query rectangle takes $O(\log^3 n + k)$ time and $O(n\log^3 n)$ space, where $n$ is the total number of polygons and $k$ is the number of polygons reported.
- Finding all polygons intersecting an orthogonal query rectangle takes $O(\log^2 m + i)$ time and $O(m\log m)$ space, where $m$ is the total number of polygon segments and $i$ is the number of segments of the polygons reported.
Now, I feel that I can finally consider the z-order. I will assume we want to query for shapes fully contained in our rectangle, since the reasoning and conclusion is similar for the other case.
I do not know why you need the objects sorted on z-order, but I assume you want to iterate over all reported shapes in their z-order (If this is not what you want, there might be a far better method than sorting)
If we query as discussed and sort the found rectangles, we get a running time of $O(\log^3 n + k\log k)$. Generally, $k$ is small compared to $n$, so the additional log factor does not do that much so you might as well sort. In fact, $n\log n$ is usually the best you can get for an geometric algorithm that is not a variant of point location/range search.
If $k$ is not small compared to $n$, you're probably better of not running all this complicated range searching and just run the trivial $O(n)$-algorithm that traverses all shapes in their z-order, while checking in constant time if it is in your query range. (I ignore sorting all shapes, as you only need to do that once)
If range searching is your problem, sorting probably isn't and if sorting is your problem, range searching probably isn't.
Although there are probably cases where they both are problem (e.g. zooming in maps), getting a good solution for this case could be rather tricky. I hope that if you learned anything from this answer, it is that range searching can already be quite a tricky business, so I'll leave it at this for now. It might be interesting to look if this specific topic has some useful literature.