From a57b3a674948ea3239925258e167a259e7a1e55b Mon Sep 17 00:00:00 2001
From: apryet <alexandre.pryet@ensegid.fr>
Date: Thu, 4 Apr 2024 14:19:30 -0400
Subject: [PATCH] Specifying raster file to zonal_stats in
 setup_basic_stress_data() lead to proj issues. Bypassed when source raster is
 loaded with rasterio and affine transform specified

---
 mfsetup/bcs.py | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/mfsetup/bcs.py b/mfsetup/bcs.py
index 7869c2a8..ba2b1983 100644
--- a/mfsetup/bcs.py
+++ b/mfsetup/bcs.py
@@ -166,7 +166,11 @@ def setup_basic_stress_data(model, shapefile=None, csvfile=None,
                 if meta['transform'][0] > m.modelgrid.delr[0]:
                     all_touched = True
                 stat = entry['stat']
-                results = zonal_stats(polygons, filename, stats=stat,
+                # load raster and specify affine transform to avoid issues with proj
+                with rasterio.open(filename) as src:
+                    affine = src.transform
+                    array = src.read(1)
+                results = zonal_stats(polygons, array, affine=affine, stats=stat,
                                     all_touched=all_touched)
                 #values = np.ones((m.nrow * m.ncol), dtype=float) * np.nan
                 #values[cells_with_bc] = np.array([r[stat] for r in results])