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": [