Fix raster merge.

spaceml-org · Dec 20, 2023 · e95c0d5 · e95c0d5
1 parent 87f7e00
commit e95c0d5
Show file tree

Hide file tree

Showing 8 changed files with 1,304 additions and 129 deletions.
diff --git a/devel_merge/Boulia_Patches.gdf b/devel_merge/Boulia_Patches.gdf
diff --git a/devel_merge/Grid_Spike.ipynb b/devel_merge/Grid_Spike.ipynb
@@ -1,5 +1,13 @@
 {
  "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "d1ad0da9",
+   "metadata": {},
+   "source": [
+    "# Creating a raster template from a patch shape"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -187,7 +195,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "d138dfc4",
+   "id": "20f65fcb",
    "metadata": {},
    "source": [
     "## Select the patch and grid"
@@ -229,6 +237,7 @@
    "outputs": [],
    "source": [
     "# Select the geometry of the patch\n",
+    "# https://epsg.io/4326\n",
     "query = (f\"SELECT patch_name, ST_AsText(geometry) \"\n",
     "         f\"FROM grid_loc \"\n",
     "         f\"WHERE patch_name = %s;\")\n",
@@ -246,22 +255,44 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "df866d40",
+   "id": "e0583fdd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def get_transform_from_geom(geom, num_pixels=2500):\n",
+    "    \"\"\"\n",
+    "    Extract an Affine transformation from a geometry, given a pixel scale\n",
+    "    \"\"\"\n",
+    "    \n",
+    "    west = geom.bounds[0]\n",
+    "    south = geom.bounds[1]\n",
+    "    east = geom.bounds[2]\n",
+    "    north = geom.bounds[3]\n",
+    "    xsize = abs(east - west)/num_pixels\n",
+    "    ysize = abs(north - south)/num_pixels\n",
+    "\n",
+    "    # Get the transformation from_origin(west, north, xsize, ysize)\n",
+    "    transform = from_origin(west, north, xsize, ysize)\n",
+    "    return transform"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "20e10fb0",
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Create a transformation given the pixel scale and coordinates\n",
+    "# Calculate the transformation matrix\n",
+    "num_pixels = 2500\n",
     "geom = grid_gdf.loc[0, 'geometry']\n",
-    "npix = 2500\n",
-    "west = geom.bounds[0]\n",
-    "south = geom.bounds[1]\n",
-    "east = geom.bounds[2]\n",
-    "north = geom.bounds[3]\n",
-    "xsize = abs(east - west)/npix\n",
-    "ysize = abs(north - south)/npix\n",
+    "transform = get_transform_from_geom(geom, num_pixels)\n",
     "\n",
-    "# Get the transformation from_origin(west, north, xsize, ysize)\n",
-    "transform = from_origin(west, north, xsize, ysize)"
+    "# At this point we have everything we need for the reprojection\n",
+    "#dst_transform\n",
+    "#dst_height = num_pixels\n",
+    "#dst_width = num_pixels\n",
+    "#dst_crs = CRS.from_epsg(\"4326\")"
    ]
   },
   {
@@ -273,9 +304,9 @@
    "source": [
     "# Create a template raster file\n",
     "#https://gis.stackexchange.com/questions/279953/numpy-array-to-gtiff-using-rasterio-without-source-raster\n",
-    "arr = np.random.randint(5, size=(npix, npix)).astype(np.uint16)\n",
-    "\n",
-    "new_dataset = rasterio.open('GRID22412_template.tif', 'w', \n",
+    "arr = np.random.randint(5, size=(num_pixels, num_pixels)).astype(np.uint16)\n",
+    "filename = 'GRID22412_template.tif'\n",
+    "new_dataset = rasterio.open(filename, 'w', \n",
     "                            driver='GTiff',\n",
     "                            height = arr.shape[0], \n",
     "                            width = arr.shape[1],\n",
@@ -290,7 +321,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "931e4e90",
+   "id": "807258af",
    "metadata": {},
    "outputs": [],
    "source": []

diff --git a/devel_merge/Merge_Spike.ipynb b/devel_merge/Merge_Spike.ipynb
@@ -28,6 +28,7 @@
     "from shapely.geometry import box, Point, MultiPoint\n",
     "from shapely.ops import nearest_points\n",
     "\n",
+    "from typing import Tuple, Callable, Any, Optional\n",
     "\n",
     "from ml4floods.data.worldfloods import dataset\n",
     "from ml4floods.data import utils\n",
@@ -127,7 +128,8 @@
    "outputs": [],
    "source": [
     "# EDIT THIS CELL: All work is conducted under a unique session name\n",
-    "session_name = \"boulia_test\""
+    "session_name = \"boulia_test\"\n",
+    "session_name = \"EMSR586\""
    ]
   },
   {
@@ -161,10 +163,10 @@
     "flood_start_date = datetime.combine(flood_start_date, midnight).replace(tzinfo=tz)\n",
     "flood_end_date = session_df.iloc[0][\"flood_date_end\"]\n",
     "flood_end_date = datetime.combine(flood_end_date, midnight).replace(tzinfo=tz)\n",
-    "ref_start_date = session_df.iloc[0][\"ref_date_start\"]\n",
-    "ref_start_date = datetime.combine(ref_start_date, midnight).replace(tzinfo=tz)\n",
-    "ref_end_date = session_df.iloc[0][\"ref_date_end\"]\n",
-    "ref_end_date = datetime.combine(ref_end_date, midnight).replace(tzinfo=tz)\n",
+    "#ref_start_date = session_df.iloc[0][\"ref_date_start\"]\n",
+    "#ref_start_date = datetime.combine(ref_start_date, midnight).replace(tzinfo=tz)\n",
+    "#ref_end_date = session_df.iloc[0][\"ref_date_end\"]\n",
+    "#ref_end_date = datetime.combine(ref_end_date, midnight).replace(tzinfo=tz)\n",
     "bucket_uri = session_df.iloc[0][\"bucket_uri\"]\n",
     "\n",
     "# Query the selected grid positions and LGAs\n",
@@ -207,10 +209,15 @@
     "\n",
     "We will use the following patches during development and testing:\n",
     "\n",
+    "**Boulia**\n",
     "* GRID22075\n",
     "* GRID22244\n",
     "* GRID21909\n",
-    "* GRID21910"
+    "* GRID21910\n",
+    "\n",
+    "**EMSR586**\n",
+    "* GRID32675\n",
+    "* GRID32844"
    ]
   },
   {
@@ -221,7 +228,9 @@
    "outputs": [],
    "source": [
     "# Set the pactch to be displayed\n",
-    "patch_name = \"GRID22075\"\n",
+    "#patch_name = \"GRID22075\"\n",
+    "patch_name = \"GRID32675\"\n",
+    "\n",
     "\n",
     "# Select the list of flood-maps during the time range\n",
     "# Order by satellite so as S2 are first\n",
@@ -346,16 +355,6 @@
     "        src.close()"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "f20276b1",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "reproj_file_list"
-   ]
-  },
   {
    "cell_type": "markdown",
    "id": "cb66afa3",
@@ -426,6 +425,14 @@
     "plt.show()"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c10230d3",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
   {
    "cell_type": "code",
    "execution_count": null,