Widen read_vrt buffer to fit all selected band dtypes (#1701)

* Widen read_vrt output buffer to fit all selected band dtypes

read_vrt allocated the multi-band output buffer from selected_bands[0].dtype
only. Each band's source array was then assigned with
result[..., band_idx] = src_arr[...], which silently casts the source to
the narrow buffer dtype. A Float32 band 1 after a Byte band 0 returned
uint8 with the float values truncated; a Byte band with
<ScaleRatio>0.5</ScaleRatio> returned uint8 with the post-scale
fractional part lost.

Compute the effective per-band dtype (declared dtype, or float64 when
any source has scale or offset, matching the existing promotion at
_vrt.py L562-565) and take np.result_type across all selected bands
before allocating the buffer. The single-band branch follows the same
logic so a single-band scaled VRT also widens. All-integer VRTs without
scaling stay integer, so memory is not blown up for the common case.

Fixes #1696

* Address Copilot review on PR #1701: drop line-number cites, guard empty VRT

Three issues raised in review:

* `_vrt.py` allocation comment cited specific line numbers ("see L562-565")
  for the ComplexSource scaling block. Line numbers drift; replace with a
  named reference to the `# Apply ComplexSource scaling` block.
* The same comment claimed mixes "widen to float64". `np.result_type` may
  also produce `float32` (e.g. NumPy 2.x on `uint8 + float32`) or
  `complex128` when complex bands are present. Reword to describe the
  common-dtype rule and list the typical outcomes.
* `test_vrt_multiband_dtype_1696.py` module docstring and one test
  docstring cited `_vrt.py` L327-334 / L562-565. Replace with named
  references (`parse_vrt` ComplexSource branch, `# Apply ComplexSource
  scaling` block) that survive future edits.
* A VRT with zero `<VRTRasterBand>` elements made `np.result_type(*[])`
  raise the generic "at least one array or dtype is required" error. Add
  an explicit `if not selected_bands` guard that raises a clear
  ValueError, plus a regression test asserting the new message.

Author: brendancol

xrspatial/geotiff/_vrt.py

@@ -522,15 +522,54 @@ def read_vrt(vrt_path: str, *, window=None,
     else:
         selected_bands = vrt.bands
 
-    # Allocate output
+    # Allocate output.
+    #
+    # The output buffer dtype must be wide enough to hold every selected
+    # band losslessly. A naive ``selected_bands[0].dtype`` would silently
+    # truncate any wider band's values when the per-band placement at the
+    # ``result[...] = src_arr[...]`` line below casts to the buffer dtype.
+    # Two sources of widening matter:
+    #
+    # * Heterogeneous declared dtypes across bands (e.g. ``Byte`` + ``Float32``).
+    # * ``ComplexSource`` ``ScaleRatio`` / ``ScaleOffset`` promote source
+    #   values to ``float64`` before placement (see the ``# Apply
+    #   ComplexSource scaling`` block later in this function); the
+    #   destination has to be float-typed too, otherwise the fractional
+    #   part is lost.
+    #
+    # ``np.result_type`` produces the common dtype that holds every
+    # contributing dtype without loss. For an all-integer VRT it stays
+    # integer; for mixes it widens to the common dtype (typically
+    # ``float64`` when integer and floating-point bands mix, but it can
+    # be a narrower float -- e.g. ``float32`` under NumPy 2.x for
+    # ``uint8`` + ``float32`` -- or ``complex128`` when complex bands
+    # are present). See issue #1696.
+    #
+    # Guard against a malformed VRT with zero ``<VRTRasterBand>``
+    # elements: ``np.result_type()`` with no args raises a generic
+    # "at least one array or dtype is required" message that gives the
+    # caller no hint about the underlying cause.
+    if not selected_bands:
+        raise ValueError(
+            "VRT has no <VRTRasterBand> elements; cannot determine "
+            "output dtype"
+        )
+    effective_dtypes = []
+    for vrt_band in selected_bands:
+        eff = vrt_band.dtype
+        for src in vrt_band.sources:
+            scaled = src.scale is not None and src.scale != 1.0
+            offset = src.offset is not None and src.offset != 0.0
+            if scaled or offset:
+                eff = np.dtype(np.float64)
+                break
+        effective_dtypes.append(eff)
+    dtype = np.result_type(*effective_dtypes)
+    fill = np.nan if dtype.kind in ('f', 'c') else 0
     if len(selected_bands) == 1:
-        dtype = selected_bands[0].dtype
-        result = np.full((out_h, out_w), np.nan if dtype.kind == 'f' else 0,
-                         dtype=dtype)
+        result = np.full((out_h, out_w), fill, dtype=dtype)
     else:
-        dtype = selected_bands[0].dtype
-        result = np.full((out_h, out_w, len(selected_bands)),
-                         np.nan if dtype.kind == 'f' else 0, dtype=dtype)
+        result = np.full((out_h, out_w, len(selected_bands)), fill, dtype=dtype)
 
     for band_idx, vrt_band in enumerate(selected_bands):
         nodata = vrt_band.nodata