to_geotiff(pack=True) re-emits the source's full per-band SCALE/OFFSET after a band-subset read; round trip corrupts values

[brendancol]

Found during a release-readiness audit of the geotiff package.

_pack() deliberately keeps attrs['gdal_metadata'] so the packed file still declares how to unpack (_attrs.py:1726-1729, it only pops scale_factor/add_offset/mask_and_scale_dtype). That is correct for full-file round trips, but after a band-subset read the kept metadata describes bands the output file no longer has.

Repro with a 2-band uint16 file carrying per-band SCALE of 0.1 (band 0) and 0.2 (band 1), raw value 50 in band 1:

da = open_geotiff("two_band.tif", band=1, unpack=True)   # 50 * 0.2 = 10.0, attrs scale_factor=0.2 - correct
to_geotiff(da, "packed.tif", pack=True)                   # reverses with 0.2, writes raw 50 - correct so far
open_geotiff("packed.tif", unpack=True)                   # raises MixedBandMetadataError
open_geotiff("packed.tif", unpack=True, band=0)           # applies band 0's stale SCALE=0.1: returns 5.0, not 10.0

The single-band output file carries both ('SCALE',0)=0.1 and ('SCALE',1)=0.2. Reading it back with unpack=True raises MixedBandMetadataError, and band=0 (the only valid index) applies the wrong scale with no error. Both outcomes contradict the pack docstring, which promises the re-packed file "unpacks to the same values".

Suggested fix: when _pack runs on a band-subset read, rewrite the per-band GDAL_METADATA down to the bands actually present (re-indexed from 0), or refuse the band-subset + pack combination with a clear error.