crop, default settings: Changes crop settings to maximum size allowed…

… found by manual work in colab. Corrects out of memory error with netcdf rotations with temporary solution to an initial crop before grid correction.
NOAA-HES-Capstone · Mar 18, 2024 · 405fdc8 · 405fdc8
1 parent ec3683d
commit 405fdc8
Show file tree

Hide file tree

Showing 4 changed files with 60 additions and 28 deletions.
diff --git a/config/preprocess_netcdf.toml b/config/preprocess_netcdf.toml
@@ -4,10 +4,10 @@ input_dir = "ims_1km"
 
 # Change this output_dir name to reflect the window size and center coordinates selected
 # Example: _50000x_15000y
-output_dir = "ims_netcdf_1km_cropped_1_000km_window"
+output_dir = "ims_netcdf_1km_cropped_2_000km_window"
 
-# Final shape will be (window_size*2, window_size*2)
-window_size = 1000 #km
+# Final shape will be (window_size, window_size)
+window_size = 2000 #km
 
 # Longitude, latitude coordinates IN DEGREES for center of cropped window. 
 # Defaults to center of current grid system if either is "None"

diff --git a/icedyno/preprocess/crop.py b/icedyno/preprocess/crop.py
@@ -61,29 +61,40 @@ def run(self) -> None:
             if os.path.exists(output_filename):
                 print(cdf_filepath, "already on disk, skipping...")
 
-            ds = xr.open_dataset(cdf_filepath, engine="h5netcdf")
+            with xr.open_dataset(cdf_filepath, engine="h5netcdf") as ds:
+                # If specified, center the window on the lat/lon coordinates provided. Otherwise, center on middle of grid.
+                if (self.center_latitude != "None") and (
+                    self.center_longitude != "None"
+                ):
+                    # Project the lat/lon coordinates to the polar stereographic coordinates.
+                    x, y = icedyno.preprocess.geolocation.polar_lonlat_to_xy(
+                        longitude=self.center_longitude, latitude=self.center_latitude
+                    )
+                else:
+                    x = np.min(np.abs(ds.x))
+                    y = np.min(np.abs(ds.y))
+
+                window = self.window_size * 1000  # from km to meters
+
+                ds = ds.sel(
+                    x=slice(x - window * 5, x + window * 5),
+                    y=slice(y - window * 5, y + window * 5),
+                )
 
-            # Correct the netCDF files with a 90 degree rotation so the xarray x,y grid matches polar stereographic
-            ds = rotate_netcdf(ds)
+                # Correct the netCDF files with a 90 degree rotation so the xarray x,y grid matches polar stereographic
+                ds = rotate_netcdf(ds)
 
-            # If specified, center the window on the lat/lon coordinates provided. Otherwise, center on middle of grid.
-            if (self.center_latitude != "None") and (self.center_longitude != "None"):
-                # Project the lat/lon coordinates to the polar stereographic coordinates.
-                x, y = icedyno.preprocess.geolocation.polar_lonlat_to_xy(
-                    longitude=self.center_longitude, latitude=self.center_latitude
+                cropped_ds = ds.sel(
+                    x=slice(x - window // 2, x + window // 2),
+                    y=slice(y - window // 2, y + window // 2),
+                )
+                assert np.allclose(
+                    cropped_ds.IMS_Surface_Values.values.shape,
+                    (1, self.window_size, self.window_size),
                 )
-            else:
-                x = np.min(np.abs(ds.x))
-                y = np.min(np.abs(ds.y))
-
-            window = self.window_size * 1000  # from km to meters
 
-            cropped_ds = ds.sel(
-                x=slice(x - window // 2, x + window // 2),
-                y=slice(y - window // 2, y + window // 2),
-            )
-            # Write the cropped data to a new NetCDF file
-            cropped_ds.to_netcdf(output_filename, engine="h5netcdf")
+                # Write the cropped data to a new NetCDF file
+                cropped_ds.to_netcdf(output_filename, engine="h5netcdf")
 
 
 def rotate_netcdf(ds: xr.Dataset) -> xr.Dataset:

diff --git a/icedyno/preprocess/geolocation.py b/icedyno/preprocess/geolocation.py
@@ -139,6 +139,6 @@ def find_closest_index_in_grid(target: float) -> int:
     grid_resolution = 1000  # meters
 
     # Calculate the index of the closest number
-    index = int((target - start) // grid_resolution)
+    index = int(np.round((target - start) // grid_resolution))
 
     return index
diff --git a/workspaces/julieanna/ims_latlong.ipynb b/workspaces/julieanna/ims_latlong.ipynb
@@ -208,7 +208,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 36,
+   "execution_count": 7,
    "id": "a7115c89-efee-409f-a39a-58aef34fcc54",
    "metadata": {},
    "outputs": [],
@@ -220,7 +220,28 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 37,
+   "execution_count": 11,
+   "id": "fa6a26e5-36d8-41c8-9d8d-92e5ccfa9ca9",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "(-1078660.2392965402, 1285497.2153701815)"
+      ]
+     },
+     "execution_count": 11,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "x, y"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
    "id": "2f9998be-7484-4f8a-bca5-2291248dc6e1",
    "metadata": {},
    "outputs": [
@@ -230,7 +251,7 @@
        "(-140.0, 73.99999999912433)"
       ]
      },
-     "execution_count": 37,
+     "execution_count": 8,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -241,7 +262,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 38,
+   "execution_count": 9,
    "id": "d46e708a-355b-4374-a304-4bc8928d14bb",
    "metadata": {},
    "outputs": [],
@@ -275,7 +296,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 39,
+   "execution_count": 10,
    "id": "d736caf8-998d-47d5-8cb5-db8246e1e1ce",
    "metadata": {},
    "outputs": [