-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
10 changed files
with
3,396 additions
and
3,341 deletions.
There are no files selected for viewing
Empty file.
Empty file.
Empty file.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,10 @@ | ||
|
|
||
| #!/bin/bash | ||
|
|
||
| mkdir -p calfin | ||
| wget -nc https://datadryad.org/stash/downloads/file_stream/458630 | ||
| mv 458630 calfin/ | ||
| cd calfin | ||
| unzip 458630 | ||
| rm 458630 | ||
| cd .. |
Binary file not shown.
Large diffs are not rendered by default.
Oops, something went wrong.
Large diffs are not rendered by default.
Oops, something went wrong.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,281 @@ | ||
| { | ||
| "cells": [ | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "6ba2c915-a0cf-4672-ad58-839841728c22", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "# Create a retreat mask from basins and Calfin" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "91ea63ac-a54a-4608-a4fa-0046e9bb40d4", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "Make sure all data sets have been downloaded prior to run this notetook" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "5bdf99df-94cc-40b4-af2a-51fcefb6aca8", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "import geopandas as gp\n", | ||
| "gp.options.io_engine = \"pyogrio\"\n", | ||
| "import pandas as pd\n", | ||
| "from tqdm.auto import tqdm\n", | ||
| "import xarray as xr\n", | ||
| "from geocube.api.core import make_geocube\n", | ||
| "import geocube\n", | ||
| "import numpy as np\n", | ||
| "from joblib import Parallel, delayed\n", | ||
| "from pathlib import Path\n", | ||
| "from typing import Union" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "4654d547-0358-4f48-9b16-617f10b28c14", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "## Prepare the IMBIE basins. Load, make valid geometry" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "2a0005de-da43-4563-8283-131931ec86bf", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "fronts_path = Path(\"../data/fronts\")\n", | ||
| "fronts_path.mkdir(parents=True, exist_ok=True)\n", | ||
| "intermediate_path = Path(\"../data/intermediate\")\n", | ||
| "intermediate_path.mkdir(parents=True, exist_ok=True)\n", | ||
| "basins_path = Path(\"../data/basins/\")\n", | ||
| "calfin_path = Path(\"../data/calfin\")" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "52f8bca6-a5ef-4457-99c0-ed1586cc01d8", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "imbie = gp.read_file(basins_path / Path(\"GRE_Basins_IMBIE2_v1.3.shp\")).to_crs(\"EPSG:3413\")\n", | ||
| "valid_geometry = imbie.make_valid()\n", | ||
| "imbie.geometry = valid_geometry\n", | ||
| "imbie = imbie[imbie[\"SUBREGION1\"] != \"ICE_CAP\"].dissolve()\n", | ||
| "valid_geometry = imbie.make_valid()\n", | ||
| "imbie.geometry = valid_geometry\n", | ||
| "imbie_gp = gp.GeoDataFrame(geometry=imbie.geometry, crs=\"EPSG:3413\")\n", | ||
| "imbie_gp.to_file(basins_path / Path(\"imbie.gpkg\"))" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "b5ae1805-1f7b-415d-a4ec-e3e8c4b3a17b", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "shelves = gp.read_file(basins_path / Path(\"GRE_Basins_shelf_extensions.gpkg\")).to_crs(\"EPSG:3413\").dissolve()\n", | ||
| "valid_geometry = shelves.make_valid()\n", | ||
| "shelves.geometry = valid_geometry\n", | ||
| "imbie = imbie.union(shelves)\n", | ||
| "imbie_gp = gp.GeoDataFrame(geometry=imbie.geometry, crs=\"EPSG:3413\")\n", | ||
| "imbie_gp.to_file(basins_path / Path(\"imbie_shelves.gpkg\"))" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "a54b0bb2-0030-4a5a-a0c4-f18123e3b5ca", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "## Prepare the IMBIE basins. Load, make valid, and dissolve" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "c96c8167-bb81-4fa9-bbac-36feac8045d3", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "calfin = gp.read_file(calfin_path / Path(\"termini_1972-2019_Greenland_closed_v1.0.shp\")).to_crs(\"EPSG:3413\")\n", | ||
| "valid_geometry = calfin.make_valid()\n", | ||
| "calfin.geometry = valid_geometry" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "881bda8d-fb39-46b6-92fd-19fd4c369d86", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "## Merge dissolved IMBIE with closed, dissolved Calfin polygons\n", | ||
| "\n", | ||
| "This creates the outline from which we are going to subtract front positions" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "f1d4117c-e51f-446b-9302-2e8d891e086f", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "calfin_dissolved = calfin.dissolve()\n", | ||
| "calfin_dissolved.to_file(intermediate_path / Path(\"calfin_dissolved.gpkg\"))" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "44b9712c-e7a7-47c7-a033-28fd160a88dd", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "union = imbie.union(calfin_dissolved)\n", | ||
| "union_gp = gp.GeoDataFrame(geometry=union, crs=\"EPSG:3413\").dissolve()\n", | ||
| "union_gp.to_file(intermediate_path / Path(\"union.gpkg\"))" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "d8afb4c3-8599-455f-bf7f-87d410d6f960", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "# Create the geometry / bounding box\n", | ||
| "# This was taken from\n", | ||
| "# https://gis.stackexchange.com/questions/436243/how-to-properly-use-geocube-geom-parameter-with-crs-other-than-4326\n", | ||
| "bbox = (-653000.0, -3384350., 879700.0, -632750.0)\n", | ||
| "geom = geocube.geo_utils.geobox.mapping(geocube.geo_utils.geobox.box(*bbox))\n", | ||
| "geom[\"crs\"] = {\"properties\": {\"name\": \"EPSG:3413\"}}" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "7d237f74-0c4e-42d1-848d-07c7ec5c0432", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "def create_ds(ds1, ds2, date, geom, crs: str = \"EPSG:3413\", grid: float = 150, out_path: Path = Path(\"../data/fronts\")):\n", | ||
| " \n", | ||
| " front = ds1.dissolve().buffer(5)\n", | ||
| " front.to_file(out_path / Path(f\"g{grid}m_terminus_{date.year}-{date.month}.gpkg\"))\n", | ||
| " diff = ds2.difference(front)\n", | ||
| " n = len(diff)\n", | ||
| " diff_df = {\"land_ice_area_fraction_retreat\": np.ones(n)}\n", | ||
| " diff_gp = gp.GeoDataFrame(data=diff_df, geometry=diff, crs=crs)\n", | ||
| " diff_gp.to_file(out_path / Path(f\"g{grid}m_frontretreat_{date.year}-{date.month}.gpkg\"))\n", | ||
| " ds = make_geocube(vector_data=diff_gp, geom=geom, resolution=[grid, grid]).fillna(0)\n", | ||
| " ds = ds.expand_dims(\"time\")\n", | ||
| " ds[\"time\"] = (\"time\", [pd.to_datetime(f\"{date.year}-{date.month}-01\")], \n", | ||
| " {\n", | ||
| " \"_FillValue\": False,\n", | ||
| " \"units\": \"days since 1972-01-01\",\n", | ||
| " \"calendar\": \"gregorian_proleptic\",\n", | ||
| " \"axis\": \"T\",\n", | ||
| " \"long_name\": \"time\",\n", | ||
| " \"unlimited\": True\n", | ||
| " },\n", | ||
| " )\n", | ||
| " ds[\"time_bounds\"] = ([\"time\", \"bounds\"], [[pd.to_datetime(f\"{date.year}-{date.month}-01\") , \n", | ||
| " pd.to_datetime(f\"{date.year}-{date.month+1}-01\")]], \n", | ||
| " { \"_FillValue\": False, \"coordinates\": False\n", | ||
| " },\n", | ||
| " )\n", | ||
| " ds[\"time\"].attrs = {\"bounds\": \"time_bounds\"}\n", | ||
| " ds[\"land_ice_area_fraction_retreat\"].attrs = {\"units\": \"1\"}\n", | ||
| " fn = out_path / Path(f\"g{grid}m_frontretreat_{date.year}-{date.month}.nc\")\n", | ||
| " comp = dict(zlib=True, complevel=2)\n", | ||
| " encoding = {var: comp for var in ds.data_vars if var != \"time_bounds\"}\n", | ||
| " ds.time.encoding = {\"units\": \"days since 1972-01-01\"}\n", | ||
| " del ds.time_bounds.attrs[\"coordinates\"]\n", | ||
| " ds.time_bounds.encoding = {\"units\": \"days since 1972-01-01\", \"coordinates\": None}\n", | ||
| " ds.to_netcdf(fn, encoding=encoding)\n", | ||
| " del ds\n", | ||
| " return fn" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "a7fbc0b0-dede-423c-91ab-151a63ff78b2", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "grid = 150" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "9379554a-3ac6-4401-8758-006eae1621e8", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "calfin[\"Date\"] = pd.to_datetime(calfin[\"Date\"])\n", | ||
| "calfin_ds = calfin.set_index(\"Date\").groupby(by=pd.Grouper(freq='ME'))\n", | ||
| "nt = len(calfin_ds)\n", | ||
| "\n", | ||
| "crs: str = \"EPSG:3413\"\n", | ||
| "\n", | ||
| "result = []\n", | ||
| "i = 0\n", | ||
| "for date, ds in tqdm(calfin_ds):\n", | ||
| " if len(ds) > 0:\n", | ||
| " if i > 0:\n", | ||
| " new_geom = old_ds.dissolve().union(ds.dissolve())\n", | ||
| " ds = gp.GeoDataFrame(geometry=new_geom, crs=crs)\n", | ||
| " fn = create_ds(ds, imbie, date, geom, grid=grid, out_path=fronts_path)\n", | ||
| " result.append(fn)\n", | ||
| " old_ds = ds\n", | ||
| " i += 1\n" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "74d4e6ed-0166-4432-9130-71b10bf66e54", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "ds = xr.open_mfdataset(result, parallel=True,)\n", | ||
| "comp = dict(zlib=True, complevel=2)\n", | ||
| "encoding = {var: comp for var in ds.data_vars if var != \"time_bounds\"}\n", | ||
| "ds.time.encoding.update({\"_FillValue\": None})\n", | ||
| "ds.time_bounds.encoding.update({\"_FillValue\": None, \"coordinates\": None})\n", | ||
| "ds.to_netcdf(fronts_path / Path(f\"pism_g{grid}m_frontretreat_calfin_1972_2019.nc\"), encoding=encoding, unlimited_dims='time')" | ||
| ] | ||
| } | ||
| ], | ||
| "metadata": { | ||
| "kernelspec": { | ||
| "display_name": "Python 3 (ipykernel)", | ||
| "language": "python", | ||
| "name": "python3" | ||
| }, | ||
| "language_info": { | ||
| "codemirror_mode": { | ||
| "name": "ipython", | ||
| "version": 3 | ||
| }, | ||
| "file_extension": ".py", | ||
| "mimetype": "text/x-python", | ||
| "name": "python", | ||
| "nbconvert_exporter": "python", | ||
| "pygments_lexer": "ipython3", | ||
| "version": "3.11.6" | ||
| } | ||
| }, | ||
| "nbformat": 4, | ||
| "nbformat_minor": 5 | ||
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -21,3 +21,12 @@ pyDOE2 | |
| cftime | ||
| nc-time-axis | ||
| netCDF4 | ||
| seaborn | ||
| pyarrow | ||
| rioxarray | ||
| geocube | ||
| geopandas | ||
| cf-xarray | ||
| nco | ||
| cdo | ||
| gsl | ||
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters