Fixes #899, #902: fix dask zonal.stats() bug, add dask+cupy backend, edge-case tests (#911)

* Fixes #899: fix boolean short-circuit bug in dask zonal.stats()

The conditions `if 'mean' or 'std' or 'var' in stats_funcs` always
evaluated to True because the string 'mean' is truthy. This caused
compute_sum, compute_count, and compute_sum_squares to always be set,
wasting work on every dask zonal.stats() call regardless of which
stats were requested.

Fix: use `any(s in stats_funcs for s in (...))` for correct membership
testing. Add regression tests covering 7 stat subsets on both numpy
and dask backends to exercise each compute flag independently.

* Fixes #900: rewrite zonal.apply() with multi-backend dispatch, return semantics

Replaces the numpy-only, input-mutating apply() with a proper
multi-backend implementation (numpy, cupy, dask+numpy, dask+cupy) that
returns a new DataArray instead of mutating the input. Uses
ArrayTypeFunctionMapping for dispatch and .data.dtype for validation
to avoid materializing dask/cupy arrays.

* Fixes #902: add dask+cupy backend for zonal.stats(), add edge-case tests

Add _stats_dask_cupy() that converts dask+cupy blocks to numpy via
map_blocks(x.get()) then delegates to the existing _stats_dask_numpy
pipeline. Wire it into the ArrayTypeFunctionMapping dispatcher.

Add five new edge-case test groups (18 test cases across backends):
- all-NaN zone: documents per-backend empty-zone behavior
- single-cell zones: std/var must be 0, not NaN
- negative zone IDs: exercises sort-and-stride with negatives
- nodata wipes zone: all finite values match nodata_values
- zone in subset of blocks: zone present in only some dask chunks

Author: brendancol