Fix a bug in GeoArray.get_mapPos() which leads to too small output in case...

Fix a bug in GeoArray.get_mapPos() which leads to too small output in case mapBounds_prj is unequal to GeoArray.prj.
Signed-off-by: Daniel Scheffler <danschef@gfz-potsdam.de>

HISTORY.rst

@@ -2,6 +2,14 @@
 History
 =======
 
+0.15.1 (22.10.2021)
+-------------------
+
+* Slight improvements for GeoArray.show_map().
+* Fixed a bug in GeoArray.get_mapPos() which leads to too small output
+  in case mapBounds_prj is unequal to GeoArray.prj (!29).
+
+
 0.15.0 (27.09.2021)
 -------------------
 

geoarray/subsetting.py

@@ -228,9 +228,11 @@ def get_array_at_mapPos(arr, arr_gt, arr_prj, out_prj, mapBounds, mapBounds_prj=
         # mapBounds are expected to have the same projection like the input array
         if not prj_equal(arr_prj, mapBounds_prj):
             xmin, ymin, xmax, ymax = mapBounds
-            (xmin, ymin), (xmax, ymax) = \
-                [transform_any_prj(mapBounds_prj, arr_prj, X, Y) for X, Y in [(xmin, ymin), (xmax, ymax)]]
-            mapBounds = xmin, ymin, xmax, ymax
+            ULxy, URxy, LRxy, LLxy = \
+                [transform_any_prj(mapBounds_prj, arr_prj, X, Y)
+                 for X, Y in [(xmin, ymax), (xmax, ymax), (xmax, ymin), (xmin, ymin)]]
+            xvals, yvals = zip(ULxy, URxy, LRxy, LLxy)
+            mapBounds = min(xvals), min(yvals), max(xvals), max(yvals)
 
         out_prj = arr_prj
         out_arr, out_gt = _clip_array_at_mapPos(arr, mapBounds, arr_gt, band2clip=band2get, fillVal=fillVal)

tests/test_geoarray.py

@@ -615,11 +615,31 @@ class Test_GeoArray(TestCase):
         xmin, xmax, ymin, ymax = self.gA.box.boundsMap
         xgsd = self.gA.xgsd
         ygsd = self.gA.ygsd
-        arr_sub, gt, prj = self.gA.get_mapPos(mapBounds=(xmin + xgsd, ymin, xmax + xgsd, ymax - ygsd / 2),
-                                              mapBounds_prj=self.gA.prj)
+        mapBounds = (xmin + xgsd, ymin, xmax + xgsd, ymax - ygsd / 2)
+
+        arr_sub, gt, prj = self.gA.get_mapPos(mapBounds=mapBounds, mapBounds_prj=self.gA.prj)
+
+        self.assertIsInstance(arr_sub, np.ndarray)
+        self.assertIsInstance(gt, (tuple, list))
+        self.assertIsInstance(prj, str)
+
+        # -- validate that the result of get_mapPos is larger if mapBounds_prj is unequal to GeoArray.prj --
+
+        # get map bounds in UTM zone 32
+        xmin, ymin, xmax, ymax = mapBounds
+        ULxy, URxy, LRxy, LLxy = \
+            [transform_any_prj(self.gA.prj, 32632, X, Y)
+             for X, Y in [(xmin, ymax), (xmax, ymax), (xmax, ymin), (xmin, ymin)]]
+        xvals, yvals = zip(ULxy, URxy, LRxy, LLxy)
+        mapBounds_32 = min(xvals), min(yvals), max(xvals), max(yvals)
+
+        arr_sub_32, gt_32, prj_32 = self.gA.get_mapPos(mapBounds=mapBounds_32, mapBounds_prj=32632)
+
         self.assertIsInstance(arr_sub, np.ndarray)
         self.assertIsInstance(gt, (tuple, list))
         self.assertIsInstance(prj, str)
+        self.assertTrue(arr_sub_32.size > arr_sub.size)
+
 
     def test_read_pointData__singleCoord_oneband(self):
         vals_z = self.gA.read_pointData(mapXY_points=np.array([[366015, 5916585]]),