✅ Make selective_tile unit test more tacit and fix it

To make the data_prep.selective_tile more understandable, I've 'simplified' the doctest with a diagonal array, and actually used more geographically correct (corner-based) pixel bounds. Also updated test_data_prep.py to use the v0.7.0 release grid instead of v0.4.0.

Issue with rasterio.open not having proper affine transformation on netcdf files solved via the most ludricrous method of all - importing xarray before rasterio... Will submit a bug report after this, but in the meantime, the code can stay mostly intact, phew! Next step is to possibly remove salem, and refactor the selective_tile function a bit more actually, by using xarray.open_rasterio and the xarray.DataArray.sel by xy slice method.

Author: weiji14

data_prep.ipynb

@@ -50,6 +50,10 @@
     "import yaml\n",
     "import zipfile\n",
     "\n",
+    "# need to import before rasterio\n",
+    "import xarray as xr\n",
+    "import salem\n",
+    "\n",
     "import geopandas as gpd\n",
     "import pygmt as gmt\n",
     "import IPython.display\n",
@@ -64,9 +68,6 @@
     "import shapely.geometry\n",
     "import skimage.util.shape\n",
     "import tqdm\n",
-    "import xarray as xr\n",
-    "\n",
-    "import salem\n",
     "\n",
     "print(\"Python       :\", sys.version.split(\"\\n\")[0])\n",
     "print(\"Geopandas    :\", gpd.__version__)\n",
@@ -1171,7 +1172,7 @@
     "def selective_tile(\n",
     "    filepath: str,\n",
     "    window_bounds: list,\n",
-    "    padding: int = 0,\n",
+    "    padding: int = 0,  # in projected coordinate system units\n",
     "    out_shape: tuple = None,\n",
     "    gapfill_raster_filepath: str = None,\n",
     ") -> np.ndarray:\n",
@@ -1182,21 +1183,21 @@
     "    some desired shape/resolution.\n",
     "\n",
     "    >>> xr.DataArray(\n",
-    "    ...     data=np.random.RandomState(seed=42).rand(64).reshape(8, 8),\n",
-    "    ...     coords={\"x\": np.arange(8), \"y\": np.arange(8)},\n",
-    "    ...     dims=[\"x\", \"y\"],\n",
+    "    ...     data=np.flipud(m=np.diag(v=np.arange(8))).astype(dtype=np.float32),\n",
+    "    ...     coords={\"y\": np.linspace(7, 0, 8), \"x\": np.linspace(0, 7, 8)},\n",
+    "    ...     dims=[\"y\", \"x\"],\n",
     "    ... ).to_netcdf(path=\"/tmp/tmp_st.nc\", mode=\"w\")\n",
     "    >>> selective_tile(\n",
     "    ...    filepath=\"/tmp/tmp_st.nc\",\n",
-    "    ...    window_bounds=[(1.0, 4.0, 3.0, 6.0), (2.0, 5.0, 4.0, 7.0)],\n",
+    "    ...    window_bounds=[(0.5, 0.5, 2.5, 2.5), (2.5, 1.5, 4.5, 3.5)],\n",
     "    ... )\n",
     "    Tiling: /tmp/tmp_st.nc ... done!\n",
-    "    array([[[[0.18485446, 0.96958464],\n",
-    "             [0.4951769 , 0.03438852]]],\n",
+    "    array([[[[0., 2.],\n",
+    "             [1., 0.]]],\n",
     "    <BLANKLINE>\n",
     "    <BLANKLINE>\n",
-    "           [[[0.04522729, 0.32533032],\n",
-    "             [0.96958464, 0.77513283]]]], dtype=float32)\n",
+    "           [[[3., 0.],\n",
+    "             [0., 0.]]]], dtype=float32)\n",
     "    >>> os.remove(\"/tmp/tmp_st.nc\")\n",
     "    \"\"\"\n",
     "    array_list = []\n",
@@ -1226,6 +1227,7 @@
     "            )\n",
     "            assert array.ndim == 3  # check that we have shape like (1, height, width)\n",
     "            assert array.shape[0] == 1  # channel-first (assuming only 1 channel)\n",
+    "            assert not 0 in array.shape  # ensure no empty dimensions (invalid window)\n",
     "\n",
     "            try:\n",
     "                assert not array.mask.any()  # check that there are no NAN values\n",
@@ -1256,7 +1258,7 @@
     "                        f\"WARN: Tile has missing data, try passing in gapfill_raster_filepath\"\n",
     "                    )\n",
     "\n",
-    "            # assert array.shape[0] == array.shape[1]  # check that height==width\n",
+    "            # assert array.shape[1] == array.shape[2]  # check that height==width\n",
     "            array_list.append(array.data.astype(dtype=np.float32))\n",
     "        print(\"done!\")\n",
     "\n",

data_prep.py

@@ -34,6 +34,10 @@
 import yaml
 import zipfile
 
+# need to import before rasterio
+import xarray as xr
+import salem
+
 import geopandas as gpd
 import pygmt as gmt
 import IPython.display
@@ -48,9 +52,6 @@
 import shapely.geometry
 import skimage.util.shape
 import tqdm
-import xarray as xr
-
-import salem
 
 print("Python       :", sys.version.split("\n")[0])
 print("Geopandas    :", gpd.__version__)
@@ -581,7 +582,7 @@ def get_window_bounds(
 def selective_tile(
     filepath: str,
     window_bounds: list,
-    padding: int = 0,
+    padding: int = 0,  # in projected coordinate system units
     out_shape: tuple = None,
     gapfill_raster_filepath: str = None,
 ) -> np.ndarray:
@@ -592,21 +593,21 @@ def selective_tile(
     some desired shape/resolution.
 
     >>> xr.DataArray(
-    ...     data=np.random.RandomState(seed=42).rand(64).reshape(8, 8),
-    ...     coords={"x": np.arange(8), "y": np.arange(8)},
-    ...     dims=["x", "y"],
+    ...     data=np.flipud(m=np.diag(v=np.arange(8))).astype(dtype=np.float32),
+    ...     coords={"y": np.linspace(7, 0, 8), "x": np.linspace(0, 7, 8)},
+    ...     dims=["y", "x"],
     ... ).to_netcdf(path="/tmp/tmp_st.nc", mode="w")
     >>> selective_tile(
     ...    filepath="/tmp/tmp_st.nc",
-    ...    window_bounds=[(1.0, 4.0, 3.0, 6.0), (2.0, 5.0, 4.0, 7.0)],
+    ...    window_bounds=[(0.5, 0.5, 2.5, 2.5), (2.5, 1.5, 4.5, 3.5)],
     ... )
     Tiling: /tmp/tmp_st.nc ... done!
-    array([[[[0.18485446, 0.96958464],
-             [0.4951769 , 0.03438852]]],
+    array([[[[0., 2.],
+             [1., 0.]]],
     <BLANKLINE>
     <BLANKLINE>
-           [[[0.04522729, 0.32533032],
-             [0.96958464, 0.77513283]]]], dtype=float32)
+           [[[3., 0.],
+             [0., 0.]]]], dtype=float32)
     >>> os.remove("/tmp/tmp_st.nc")
     """
     array_list = []
@@ -636,6 +637,7 @@ def selective_tile(
             )
             assert array.ndim == 3  # check that we have shape like (1, height, width)
             assert array.shape[0] == 1  # channel-first (assuming only 1 channel)
+            assert not 0 in array.shape  # ensure no empty dimensions (invalid window)
 
             try:
                 assert not array.mask.any()  # check that there are no NAN values
@@ -666,7 +668,7 @@ def selective_tile(
                         f"WARN: Tile has missing data, try passing in gapfill_raster_filepath"
                     )
 
-            # assert array.shape[0] == array.shape[1]  # check that height==width
+            # assert array.shape[1] == array.shape[2]  # check that height==width
             array_list.append(array.data.astype(dtype=np.float32))
         print("done!")
 

features/steps/test_data_prep.py

@@ -62,7 +62,7 @@ def get_a_raster_grid(context, dataset_type, raster_grid):
     context.raster_grid = raster_grid
     context.filepath = os.path.join(dataset_type, raster_grid)
     url = (
-        f"https://github.com/weiji14/deepbedmap/releases/download/v0.4.0/{raster_grid}"
+        f"https://github.com/weiji14/deepbedmap/releases/download/v0.7.0/{raster_grid}"
     )
     context.data_prep.download_to_path(path=context.filepath, url=url)
 

test_ipynb.ipynb

@@ -170,43 +170,26 @@
       "ok\n",
       "Trying:\n",
       "    xr.DataArray(\n",
-      "        data=np.random.RandomState(seed=42).rand(64).reshape(8, 8),\n",
-      "        coords={\"x\": np.arange(8), \"y\": np.arange(8)},\n",
-      "        dims=[\"x\", \"y\"],\n",
+      "        data=np.flipud(m=np.diag(v=np.arange(8))).astype(dtype=np.float32),\n",
+      "        coords={\"y\": np.linspace(7, 0, 8), \"x\": np.linspace(0, 7, 8)},\n",
+      "        dims=[\"y\", \"x\"],\n",
       "    ).to_netcdf(path=\"/tmp/tmp_st.nc\", mode=\"w\")\n",
       "Expecting nothing\n",
       "ok\n",
       "Trying:\n",
       "    selective_tile(\n",
       "       filepath=\"/tmp/tmp_st.nc\",\n",
-      "       window_bounds=[(1.0, 4.0, 3.0, 6.0), (2.0, 5.0, 4.0, 7.0)],\n",
+      "       window_bounds=[(0.5, 0.5, 2.5, 2.5), (2.5, 1.5, 4.5, 3.5)],\n",
       "    )\n",
       "Expecting:\n",
       "    Tiling: /tmp/tmp_st.nc ... done!\n",
-      "    array([[[[0.18485446, 0.96958464],\n",
-      "             [0.4951769 , 0.03438852]]],\n",
-      "    <BLANKLINE>\n",
-      "    <BLANKLINE>\n",
-      "           [[[0.04522729, 0.32533032],\n",
-      "             [0.96958464, 0.77513283]]]], dtype=float32)\n",
-      "**********************************************************************\n",
-      "File \"data_prep\", line 627, in data_prep.selective_tile\n",
-      "Failed example:\n",
-      "    selective_tile(\n",
-      "       filepath=\"/tmp/tmp_st.nc\",\n",
-      "       window_bounds=[(1.0, 4.0, 3.0, 6.0), (2.0, 5.0, 4.0, 7.0)],\n",
-      "    )\n",
-      "Expected:\n",
-      "    Tiling: /tmp/tmp_st.nc ... done!\n",
-      "    array([[[[0.18485446, 0.96958464],\n",
-      "             [0.4951769 , 0.03438852]]],\n",
+      "    array([[[[0., 2.],\n",
+      "             [1., 0.]]],\n",
       "    <BLANKLINE>\n",
       "    <BLANKLINE>\n",
-      "           [[[0.04522729, 0.32533032],\n",
-      "             [0.96958464, 0.77513283]]]], dtype=float32)\n",
-      "Got:\n",
-      "    Tiling: /tmp/tmp_st.nc ... done!\n",
-      "    array([], shape=(2, 1, 0, 2), dtype=float32)\n",
+      "           [[[3., 0.],\n",
+      "             [0., 0.]]]], dtype=float32)\n",
+      "ok\n",
       "Trying:\n",
       "    os.remove(\"/tmp/tmp_st.nc\")\n",
       "Expecting nothing\n",
@@ -241,19 +224,17 @@
       "2 items had no tests:\n",
       "    data_prep\n",
       "    data_prep.parse_datalist\n",
-      "6 items passed all tests:\n",
+      "7 items passed all tests:\n",
       "   6 tests in data_prep.ascii_to_xyz\n",
       "   3 tests in data_prep.check_sha256\n",
       "   3 tests in data_prep.download_to_path\n",
       "   2 tests in data_prep.get_region\n",
       "   3 tests in data_prep.get_window_bounds\n",
+      "   3 tests in data_prep.selective_tile\n",
       "   5 tests in data_prep.xyz_to_grid\n",
-      "**********************************************************************\n",
-      "1 items had failures:\n",
-      "   1 of   3 in data_prep.selective_tile\n",
       "25 tests in 9 items.\n",
-      "24 passed and 1 failed.\n",
-      "***Test Failed*** 1 failures.\n"
+      "25 passed and 0 failed.\n",
+      "Test passed.\n"
      ]
     }
    ],