{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "55654545-a837-413e-8826-310b063ff6bd",
   "metadata": {},
   "source": [
    "# Display a DEM as a 3D surface\n",
    "O. Kaufmann, 2024."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e934d22e-9987-4037-a565-22d44bd32056",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.environ[\"TRAME_DISABLE_V3_WARNING\"] = \"1\"  # disable trame v3 warning\n",
    "import pyvista as pv\n",
    "from geometron.utils.geoservices import esri_rest_server_geotiff\n",
    "from geometron.vtk import geotiff_to_grid"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a51262da-8c25-4770-9969-7196d75c61f5",
   "metadata": {},
   "source": [
    "## Get the DEM from WalOnMap (not possible because these are images?)\n",
    "https://www.imagetest.geomatics.gov.nt.ca/arcgis/rest/services/GNWT_Basemaps/Aster_DEM_Basemap/MapServer?f=pjson"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f0519054-ea15-4b63-a6da-51a15ce64d72",
   "metadata": {},
   "outputs": [],
   "source": [
    "host = 'www.imagetest.geomatics.gov.nt.ca'\n",
    "folder = 'GNWT_Basemaps'\n",
    "service = 'Aster_DEM_Basemap'\n",
    "service_type = 'MapServer'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0a3a4afd-b835-435b-8107-efd819e3b332",
   "metadata": {},
   "outputs": [],
   "source": [
    "xmin, xmax, ymin, ymax = 620400.0, 620600.0, 626700.0, 627000.0  # map extent in epsg:31370\n",
    "size_width, size_height = 768, 768\n",
    "epsg = 3857"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "45aa31d7-b15e-4bb6-b64e-6aa1a970ed26",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "params = {'bbox': f'{xmin},{ymin},{xmax},{ymax}',\n",
    "          'size': f'{size_width}, {size_height}',\n",
    "          'bboxSR':{\"wkid\" : epsg},\n",
    "          'imageSR':{\"wkid\" : epsg},\n",
    "          'format': 'image/geotiff',\n",
    "          'transparent': 'false',\n",
    "          'layerDefs': '',\n",
    "          'layers' : '',\n",
    "          'dpi' : '100',\n",
    "          'time':'',\n",
    "          'f':'image'\n",
    "         }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "2f3b3d0d-5044-42a5-ad92-7c1e5ec31ba7",
   "metadata": {},
   "outputs": [],
   "source": [
    "image_file='./data/dem.tif'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b8f93ef6-c5ff-4a2d-bf22-1b893c15ae99",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "url:\n",
      "https://www.imagetest.geomatics.gov.nt.ca/arcgis/rest/services/GNWT_Basemaps/Aster_DEM_Basemap/MapServer/export?bbox=620400.0%2C626700.0%2C620600.0%2C627000.0&size=768%2C+768&bboxSR=%7B%27wkid%27%3A+3857%7D&imageSR=%7B%27wkid%27%3A+3857%7D&format=image%2Fgeotiff&transparent=false&layerDefs=&layers=&dpi=100&time=&f=pjson\n",
      "\n",
      "url:\n",
      "https://www.imagetest.geomatics.gov.nt.ca/arcgis/rest/services/GNWT_Basemaps/Aster_DEM_Basemap/MapServer/export?bbox=620400.0%2C626700.0%2C620600.0%2C627000.0&size=768%2C+768&bboxSR=%7B%27wkid%27%3A+3857%7D&imageSR=%7B%27wkid%27%3A+3857%7D&format=image%2Fgeotiff&transparent=false&layerDefs=&layers=&dpi=100&time=&f=image\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/su530201/.local/share/virtualenvs/geometron-VhFx4lLH/lib/python3.10/site-packages/rasterio/__init__.py:304: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.\n",
      "  dataset = DatasetReader(path, driver=driver, sharing=sharing, **kwargs)\n"
     ]
    }
   ],
   "source": [
    "esri_rest_server_geotiff(image_file=image_file, host=host,folder=folder, service=service, service_type=service_type, parameter=params, verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fb4f9d5-989a-4b1f-b6b2-f43570f07384",
   "metadata": {},
   "source": [
    "## Load DEM as a grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "a08291cd-c209-4160-9218-f751dcd39395",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ERROR:root:Cannot read tile : 0,0 from file\n",
      "/home/su530201/.local/share/virtualenvs/geometron-VhFx4lLH/lib/python3.10/site-packages/pyvista/core/utilities/fileio.py:203: UserWarning: The VTK reader `vtkTIFFReader` in pyvista reader `TIFFReader('/home/su530201/PycharmProjects/geometron/docs/source/examples/geological/data/dem.tif')` raised an errorwhile reading the file.\n",
      "\t\"Cannot read tile : 0,0 from file\"\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[7], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m dem \u001b[38;5;241m=\u001b[39m \u001b[43mgeotiff_to_grid\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43m./data/dem.tif\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mwarp\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mTrue\u001b[39;49;00m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mfactor\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m5.\u001b[39;49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/PycharmProjects/geometron/geometron/vtk/image_to_vtk.py:23\u001b[0m, in \u001b[0;36mgeotiff_to_grid\u001b[0;34m(filename, warp, factor)\u001b[0m\n\u001b[1;32m     21\u001b[0m transform \u001b[38;5;241m=\u001b[39m shapely_matrix_to_vtk_matrix(rasterio_matrix_to_shapely_matrix(src\u001b[38;5;241m.\u001b[39mread_transform()), force_3d\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m)\n\u001b[1;32m     22\u001b[0m grid \u001b[38;5;241m=\u001b[39m grid\u001b[38;5;241m.\u001b[39mtransform(transform)\n\u001b[0;32m---> 23\u001b[0m no_data \u001b[38;5;241m=\u001b[39m [n \u001b[38;5;28;01mfor\u001b[39;00m n, i \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28menumerate\u001b[39m(grid\u001b[38;5;241m.\u001b[39mpoint_data[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mTiff Scalars\u001b[39m\u001b[38;5;124m'\u001b[39m]) \u001b[38;5;28;01mif\u001b[39;00m i \u001b[38;5;129;01min\u001b[39;00m src\u001b[38;5;241m.\u001b[39mnodatavals]\n\u001b[1;32m     24\u001b[0m grid\u001b[38;5;241m.\u001b[39mpoint_data[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mTiff Scalars\u001b[39m\u001b[38;5;124m'\u001b[39m][no_data] \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mnan\n\u001b[1;32m     25\u001b[0m grid\u001b[38;5;241m.\u001b[39mset_active_scalars(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mTiff Scalars\u001b[39m\u001b[38;5;124m'\u001b[39m)\n",
      "File \u001b[0;32m~/PycharmProjects/geometron/geometron/vtk/image_to_vtk.py:23\u001b[0m, in \u001b[0;36m<listcomp>\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m     21\u001b[0m transform \u001b[38;5;241m=\u001b[39m shapely_matrix_to_vtk_matrix(rasterio_matrix_to_shapely_matrix(src\u001b[38;5;241m.\u001b[39mread_transform()), force_3d\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m)\n\u001b[1;32m     22\u001b[0m grid \u001b[38;5;241m=\u001b[39m grid\u001b[38;5;241m.\u001b[39mtransform(transform)\n\u001b[0;32m---> 23\u001b[0m no_data \u001b[38;5;241m=\u001b[39m [n \u001b[38;5;28;01mfor\u001b[39;00m n, i \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28menumerate\u001b[39m(grid\u001b[38;5;241m.\u001b[39mpoint_data[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mTiff Scalars\u001b[39m\u001b[38;5;124m'\u001b[39m]) \u001b[38;5;28;01mif\u001b[39;00m \u001b[43mi\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;129;43;01min\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43msrc\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mnodatavals\u001b[49m]\n\u001b[1;32m     24\u001b[0m grid\u001b[38;5;241m.\u001b[39mpoint_data[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mTiff Scalars\u001b[39m\u001b[38;5;124m'\u001b[39m][no_data] \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mnan\n\u001b[1;32m     25\u001b[0m grid\u001b[38;5;241m.\u001b[39mset_active_scalars(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mTiff Scalars\u001b[39m\u001b[38;5;124m'\u001b[39m)\n",
      "\u001b[0;31mValueError\u001b[0m: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()"
     ]
    }
   ],
   "source": [
    "dem = geotiff_to_grid('./data/dem.tif', warp=True, factor=5.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "0a19923b-583a-4946-a65d-3f0ff6972ba4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterio as rio\n",
    "from geometron.geometries.transforms import shapely_matrix_to_vtk_matrix, rasterio_matrix_to_shapely_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "66ecac7e-6258-4ed2-bd16-a76221e75540",
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = './data/dem.tif'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "af2014d8-f0c8-4756-891b-4e7abb2fd6f1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[253, 253, 253, ..., 253, 253, 253],\n",
       "        [253, 253, 253, ..., 253, 253, 253],\n",
       "        [253, 253, 253, ..., 253, 253, 253],\n",
       "        ...,\n",
       "        [253, 253, 253, ..., 253, 253, 253],\n",
       "        [253, 253, 253, ..., 253, 253, 253],\n",
       "        [253, 253, 253, ..., 253, 253, 253]],\n",
       "\n",
       "       [[253, 253, 253, ..., 253, 253, 253],\n",
       "        [253, 253, 253, ..., 253, 253, 253],\n",
       "        [253, 253, 253, ..., 253, 253, 253],\n",
       "        ...,\n",
       "        [253, 253, 253, ..., 253, 253, 253],\n",
       "        [253, 253, 253, ..., 253, 253, 253],\n",
       "        [253, 253, 253, ..., 253, 253, 253]],\n",
       "\n",
       "       [[253, 253, 253, ..., 253, 253, 253],\n",
       "        [253, 253, 253, ..., 253, 253, 253],\n",
       "        [253, 253, 253, ..., 253, 253, 253],\n",
       "        ...,\n",
       "        [253, 253, 253, ..., 253, 253, 253],\n",
       "        [253, 253, 253, ..., 253, 253, 253],\n",
       "        [253, 253, 253, ..., 253, 253, 253]],\n",
       "\n",
       "       [[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": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "with rio.open(filename, mode='r') as src:\n",
    "    img = src.read()\n",
    "img"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "869065a2-6ca5-4bd4-bb69-38b91487450e",
   "metadata": {},
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "'numpy.ndarray' object has no attribute 'cast_to_structured_grid'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[13], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m grid \u001b[38;5;241m=\u001b[39m \u001b[43mimg\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcast_to_structured_grid\u001b[49m()\n\u001b[1;32m      2\u001b[0m src \u001b[38;5;241m=\u001b[39m rio\u001b[38;5;241m.\u001b[39mopen(filename)\n\u001b[1;32m      3\u001b[0m transform \u001b[38;5;241m=\u001b[39m shapely_matrix_to_vtk_matrix(rasterio_matrix_to_shapely_matrix(src\u001b[38;5;241m.\u001b[39mread_transform()), force_3d\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m)\n",
      "\u001b[0;31mAttributeError\u001b[0m: 'numpy.ndarray' object has no attribute 'cast_to_structured_grid'"
     ]
    }
   ],
   "source": [
    "grid = img.cast_to_structured_grid()\n",
    "src = rio.open(filename)\n",
    "transform = shapely_matrix_to_vtk_matrix(rasterio_matrix_to_shapely_matrix(src.read_transform()), force_3d=True)\n",
    "grid = grid.transform(transform)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61ebd50b-5622-45ee-9a30-2126d08833ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "no_data = [n for n, i in enumerate(grid.point_data['Tiff Scalars']) if i in src.nodatavals]\n",
    "grid.point_data['Tiff Scalars'][no_data] = np.nan\n",
    "grid.set_active_scalars('Tiff Scalars')\n",
    "if warp:\n",
    "    grid = grid.warp_by_scalar('Tiff Scalars', factor=factor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4200848d-1275-488a-bd6e-879e09fcce21",
   "metadata": {},
   "outputs": [],
   "source": [
    "grid.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bdba8b64-5a3b-4aed-ac31-50ca35e6507d",
   "metadata": {},
   "outputs": [],
   "source": [
    "grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3fe2ebca-7615-4c2a-9352-2b34f912f61f",
   "metadata": {},
   "outputs": [],
   "source": [
    "grid.set_active_scalars('Tiff Scalars')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc0e3487-1d3b-47c0-a66c-21d2223c7b6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "grid.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8dc94baa-8265-4ee8-a3df-51419b168300",
   "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}