{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "aaa66657-66d4-46c1-a43e-a4512739ebd5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#import required libraries\n",
    "import rasterio\n",
    "import matplotlib.pyplot as plt\n",
    "from rasterio.plot import show"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2c916e8a-a7d7-412a-9f5b-3d3f19ecca9e",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\saulm\\anaconda3\\Lib\\site-packages\\rasterio\\__init__.py:317: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.\n",
      "  dataset = DatasetReader(path, driver=driver, sharing=sharing, **kwargs)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<open DatasetReader name='data/Peligros_Geologicos.jpg' mode='r'>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#open ungeoreferenced raster\n",
    "unRefRaster = rasterio.open('data/Peligros_Geologicos.jpg')\n",
    "unRefRaster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "733fbd4e-1c03-4549-8757-3db0bd671dd7",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[255, 255, 255, ..., 255, 255, 255],\n",
       "       [255, 255, 255, ..., 255, 255, 255],\n",
       "       [255, 255, 255, ..., 255, 255, 255],\n",
       "       ...,\n",
       "       [255, 255, 255, ..., 255, 255, 255],\n",
       "       [255, 255, 255, ..., 255, 255, 255],\n",
       "       [255, 255, 255, ..., 255, 255, 255]], dtype=uint8)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#show raster band values\n",
    "unRefRaster.read(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "229ad982-1244-484f-9ad6-d33909b88a75",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 158636 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Axes: >"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#show raster\n",
    "show(unRefRaster)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "59993d5a-ea22-4c97-83a8-6507f0f22d3f",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4133, 2922)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#show raster shape\n",
    "unRefRaster.read(1).shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7880f183-135e-4add-b8e9-22e069af2415",
   "metadata": {},
   "source": [
    "## Insert control points from coordinates captured in paint\n",
    "\n",
    "Control point 1\n",
    "![Alt text](data/GCP1.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "693b8992-607f-44c4-a451-3486d899b01a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GroundControlPoint(row=368, col=190, x=-81, y=-1, id='c13ed452-6751-4d57-90a8-fe7de053e751')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "point1 = rasterio.control.GroundControlPoint(row=368, col=190, x=-81, y=-1)\n",
    "point1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96b2929b-739d-4b44-a221-8945e85b62d3",
   "metadata": {},
   "source": [
    "Control point 2\n",
    "![Alt text](data/GCP2.png)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "9270ee83-52d2-4997-98c8-7b50883cce43",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GroundControlPoint(row=3497, col=239, x=-81, y=-16, id='3cb08473-2784-440b-8c31-3f30403c0305')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "point2 = rasterio.control.GroundControlPoint(row=3497, col=239, x=-81, y=-16)\n",
    "point2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa6d4bd3-fed1-4cd3-9ccb-4fa22c923a58",
   "metadata": {},
   "source": [
    "Control point 3\n",
    "![Alt text](data/GCP3.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "98516fd4-1a75-40d3-ad25-da6a797d78a9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GroundControlPoint(row=3706, col=2645, x=-69, y=-17, id='45aae98b-fdce-442c-8e23-561f2589e111')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "point3 = rasterio.control.GroundControlPoint(row=3706, col=2645, x=-69, y=-17)\n",
    "point3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "979f45b8-5ade-40cb-bcb6-e5d9ee07ac8e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[GroundControlPoint(row=368, col=190, x=-81, y=-1, id='c13ed452-6751-4d57-90a8-fe7de053e751'),\n",
       " GroundControlPoint(row=3497, col=239, x=-81, y=-16, id='3cb08473-2784-440b-8c31-3f30403c0305'),\n",
       " GroundControlPoint(row=3706, col=2645, x=-69, y=-17, id='45aae98b-fdce-442c-8e23-561f2589e111')]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#list of selected gcps\n",
    "points = [point1, point2, point3]\n",
    "points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "7ea4d09f-5183-4328-a894-39923a0f36e7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Affine(0.004994325053839821, -7.821090688339848e-05, -81.92014014649648,\n",
       "       7.980704784024951e-07, -0.00479387635201452, 0.7639948641504404)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#get transformation array from points\n",
    "transformation = rasterio.transform.from_gcps(points)\n",
    "transformation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "a3155965-a0e0-4a48-b95f-a3fc79a20b50",
   "metadata": {},
   "outputs": [],
   "source": [
    "#define output raster\n",
    "outputPath = 'data/georefRaster.tif'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "235c8f90-ab1d-42cb-a52f-78a8a1570f7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#create raster and write bands\n",
    "with rasterio.open(\n",
    "    outputPath,\n",
    "    'w',\n",
    "    driver='GTiff',\n",
    "    height=unRefRaster.read(1).shape[0],\n",
    "    width=unRefRaster.read(1).shape[1],\n",
    "    count=3,\n",
    "    dtype=unRefRaster.read(1).dtype,\n",
    "    crs=rasterio.crs.CRS.from_epsg(4326),\n",
    "    transform=transformation,\n",
    ") as dst:\n",
    "    dst.write(unRefRaster.read(1), 1)\n",
    "    dst.write(unRefRaster.read(2), 2)\n",
    "    dst.write(unRefRaster.read(3), 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "f1dd4c93-0683-4954-b0ec-fe9c48402722",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 164564 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Axes: >"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#show georeferenced raster\n",
    "geoRaster = rasterio.open(outputPath)\n",
    "show(geoRaster)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1579af08-0866-4dd5-a49c-4b5be4bad47d",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
