From 405fdc826b7f09c26d9345222b2e618ca1f15f0c Mon Sep 17 00:00:00 2001 From: Julieanna Bacon Date: Mon, 18 Mar 2024 12:55:15 -0400 Subject: [PATCH] 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. --- config/preprocess_netcdf.toml | 6 ++-- icedyno/preprocess/crop.py | 49 ++++++++++++++++---------- icedyno/preprocess/geolocation.py | 2 +- workspaces/julieanna/ims_latlong.ipynb | 31 +++++++++++++--- 4 files changed, 60 insertions(+), 28 deletions(-) diff --git a/config/preprocess_netcdf.toml b/config/preprocess_netcdf.toml index cf8e343..a90f962 100644 --- 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 index 877d90b..84dae3d 100644 --- 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 index 0e5ba63..975cb15 100644 --- 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 index 22b5e9b..01f98cf 100644 --- 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": [