nsdf_geo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e3edd4eb",
   "metadata": {},
   "source": [
    "# Tutorial: Using NSDF for End-to-End Analysis of Scientific Data\n",
    "\n",
    "In this tutorial, you will learn how to:\n",
    "\n",
    "* Build a modular workflow integrating your application module with NSDF services for data visualization and analysis.\n",
    "* Load and store intermediate data from/to public storage options like Dataverse and from/to private storage such as Seal Storage.\n",
    "* Use the NSDF dashboard for data visualization and analysis. The analysis includes zooming into data for detailed analysis, cropping subregions, and saving them locally in a Python-compatible format.\n",
    "\n",
    "### **DISCLAIMER:** Please **DO NOT *Run All Cells*** because different paths can be executed leading to the same outcome"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6bf1f88d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import geotiled as gt\n",
    "from pathlib import Path\n",
    "import glob\n",
    "import os\n",
    "import shutil\n",
    "import multiprocessing\n",
    "import OpenVisus as ov\n",
    "import numpy as np\n",
    "import requests\n",
    "import json\n",
    "from matplotlib import pyplot as plt\n",
    "from tqdm import tqdm\n",
    "\n",
    "# To silence a deprecation warning.\n",
    "gt.gdal.UseExceptions()\n",
    "\n",
    "# You have have successfully prepared your environment.\n",
    "print(\"You have successfully prepared your environment.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5433cd70-319d-4e0c-b2aa-b9267fe9cc7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.exists(\"files/tif_files\"):\n",
    "    os.mkdir(\"files/tif_files\")\n",
    "\n",
    "with open(\"./files/json/dataverse.json\", \"r\") as file:\n",
    "    urls = json.load(file)\n",
    "\n",
    "\n",
    "def get_data_from_dataverse(data):\n",
    "    file_url = data.get(\"url\")\n",
    "    name_file = data.get(\"name_file\")\n",
    "    resp = requests.get(file_url)\n",
    "    with open(name_file, \"wb\") as f:\n",
    "        f.write(resp.content)\n",
    "\n",
    "\n",
    "print(\"Your download is starting...\")\n",
    "\n",
    "pool = multiprocessing.Pool(processes=len(urls))\n",
    "pool.map(get_data_from_dataverse, urls)\n",
    "\n",
    "\n",
    "print(\"You have successfully downloaded the data from Dataverse.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35d0dda7-e5ba-4be3-b6f4-865ab31bf875",
   "metadata": {},
   "outputs": [],
   "source": [
    "shapefile = [\"./files/shape_files/STATEFP_47.shp\"]  # Shapefile for Visualization\n",
    "hill = os.path.join(\"./files/tif_files/\", \"hillshading.tif\")\n",
    "pcs = os.path.join(\"./files/tif_files/\", \"elevation.tif\")\n",
    "\n",
    "pcs_array = gt.generate_img(\n",
    "    pcs,\n",
    "    downsample=5,\n",
    "    reproject_gcs=True,\n",
    "    shp_files=shapefile,\n",
    "    title=\"Elevation Data for TN @ 1 Arc-Second/30m Resolution\",\n",
    "    zunit=\"Meter\",\n",
    "    xyunit=\"Degree\",\n",
    "    ztype=\"Elevation\",\n",
    "    crop_shp=True,\n",
    ")\n",
    "hill_array = gt.generate_img(\n",
    "    hill,\n",
    "    downsample=5,\n",
    "    reproject_gcs=True,\n",
    "    shp_files=shapefile,\n",
    "    title=\"Hillshading Data for TN @ 1 Arc-Second/30m Resolution\",\n",
    "    zunit=\"Level\",\n",
    "    xyunit=\"Degree\",\n",
    "    ztype=\"Hillshading\",\n",
    "    crop_shp=True,\n",
    ")\n",
    "\n",
    "print(\"You have successfully visualized the data from Dataverse.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5adfc770",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate lat/lon min/max from tiff\n",
    "from osgeo import gdal, osr\n",
    "\n",
    "hill = os.path.join(\"./files/tif_files/\", \"hillshading.tif\")\n",
    "\n",
    "dataset = gdal.Open(hill)\n",
    "band = dataset.GetRasterBand(1)\n",
    "\n",
    "geotransform = dataset.GetGeoTransform()\n",
    "spatial_ref = osr.SpatialReference(wkt=dataset.GetProjection())\n",
    "\n",
    "target_spatial_ref = osr.SpatialReference()\n",
    "target_spatial_ref.ImportFromEPSG(4326)\n",
    "\n",
    "coord_transform = osr.CoordinateTransformation(spatial_ref, target_spatial_ref)\n",
    "ulx, xres, _, uly, _, yres = geotransform\n",
    "lrx = ulx + (dataset.RasterXSize * xres)\n",
    "lry = uly + (dataset.RasterYSize * yres)\n",
    "top_left_geo = coord_transform.TransformPoint(ulx, uly)\n",
    "bottom_right_geo = coord_transform.TransformPoint(lrx, lry)\n",
    "\n",
    "lon_1, lat_1, _ = top_left_geo\n",
    "lon_2, lat_2, _ = bottom_right_geo\n",
    "lat_min = min(lat_1, lat_2)\n",
    "lat_max = max(lat_1, lat_2)\n",
    "lon_min = min(lon_1, lon_2)\n",
    "lon_max = max(lon_1, lon_2)\n",
    "\n",
    "print(f\"Longitude range: {lat_min} to {lat_max}\")\n",
    "print(f\"Latitude range: {lon_min} to {lon_max}\")\n",
    "\n",
    "print(\"You have successfully extract geospatial metadata from the TIFFs.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "885440a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = \"idx_data/Tennessee_terrain_parameters.idx\"\n",
    "all_fields = [ov.Field(\"elevation\", \"float32\"), ov.Field(\"hillshading\", \"uint8\")]\n",
    "input_data = [np.flipud(pcs_array).copy(), np.flipud(hill_array).copy()]\n",
    "height, width = input_data[0].shape\n",
    "db = ov.CreateIdx(\n",
    "    url=filename,\n",
    "    dims=[width, height],\n",
    "    fields=all_fields,\n",
    "    arco=\"4mb\",\n",
    "    physic_box=ov.BoxNd.fromString(f\"{lat_min} {lat_max} {lon_min} {lon_max}\"),\n",
    "    time=[0, 0, \"%00000d/\"],\n",
    ")\n",
    "i = 0\n",
    "for fld in db.getFields():\n",
    "    db.write(input_data[i], field=fld)\n",
    "    i += 1\n",
    "db.compressDataset([\"zip\"])\n",
    "\n",
    "print(\"You have successfully created the IDX file with the two terrain parmaters.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5db099e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = \"Tennessee_terrain_parameters.idx\"\n",
    "HOME_DIR = \"s3://utah/nsdf/somospie/\"  # DONOT change this line\n",
    "data_dir = \"terrain_tennessee/\"\n",
    "upload_dir = HOME_DIR + data_dir\n",
    "s3_path = upload_dir.split(\"://\")[1]\n",
    "s3_path += filename\n",
    "remote_dir = (\n",
    "    \"https://maritime.sealstorage.io/api/v0/s3/\"\n",
    "    + s3_path\n",
    "    + \"?access_key=any&secret_key=any&endpoint_url=https://maritime.sealstorage.io/api/v0/s3&cached=arco\"\n",
    ")\n",
    "\n",
    "db = ov.LoadDataset(remote_dir)\n",
    "\n",
    "print(\"You have successfully downloaded the IDX files from Seal Storage.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2a340c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "read_elevation = db.read(field=\"elevation\")\n",
    "read_hillshading = db.read(field=\"hillshading\")\n",
    "read_aspect = db.read(field=\"aspect\")\n",
    "read_slope = db.read(field=\"slope\")\n",
    "\n",
    "cmap_instance = plt.get_cmap(\"inferno\")\n",
    "\n",
    "fig, axs = plt.subplots(4, 1, figsize=(10, 8))\n",
    "axs[0].set_xlim(lat_min, lat_max)\n",
    "axs[0].set_ylim(lon_min, lon_max)\n",
    "axs[0].set_title(\"Elevation - TN 30m Resolution\")\n",
    "axs[0].set_xlabel(\"Longitude (Degrees)\")\n",
    "axs[0].set_ylabel(\"Latitude (Degrees)\")\n",
    "elevation_fig = axs[0].imshow(\n",
    "    read_elevation,\n",
    "    cmap=cmap_instance,\n",
    "    origin=\"lower\",\n",
    "    vmin=30,\n",
    "    vmax=1999,\n",
    "    extent=(lat_min, lat_max, lon_min, lon_max),\n",
    ")\n",
    "\n",
    "cbar = fig.colorbar(\n",
    "    elevation_fig,\n",
    "    fraction=0.046 * read_elevation.shape[0] / read_elevation.shape[1],\n",
    "    pad=0.04,\n",
    ")\n",
    "cbar_ticks = np.linspace(30, 1999, 8)\n",
    "cbar.set_ticks(cbar_ticks)\n",
    "cbar.set_label(\"Elevation (Meters)\")\n",
    "\n",
    "axs[1].set_xlim(lat_min, lat_max)\n",
    "axs[1].set_ylim(lon_min, lon_max)\n",
    "axs[1].set_title(\"Hillshading - TN 30m Resolution\")\n",
    "axs[1].set_xlabel(\"Longitude (Degrees)\")\n",
    "axs[1].set_ylabel(\"Latitude (Degrees)\")\n",
    "hillshading_fig = axs[1].imshow(\n",
    "    read_hillshading,\n",
    "    cmap=cmap_instance,\n",
    "    origin=\"lower\",\n",
    "    vmin=0,\n",
    "    vmax=255,\n",
    "    extent=(lat_min, lat_max, lon_min, lon_max),\n",
    ")\n",
    "\n",
    "cbar = fig.colorbar(\n",
    "    hillshading_fig,\n",
    "    fraction=0.046 * read_hillshading.shape[0] / read_hillshading.shape[1],\n",
    "    pad=0.04,\n",
    ")\n",
    "cbar_ticks = np.linspace(0, 255, 8)\n",
    "cbar.set_ticks(cbar_ticks)\n",
    "cbar.set_label(\"Hillshading (Levels)\")\n",
    "\n",
    "axs[2].set_xlim(lat_min, lat_max)\n",
    "axs[3].set_ylim(lon_min, lon_max)\n",
    "axs[2].set_title(\"Aspect - TN 30m Resolution\")\n",
    "axs[2].set_xlabel(\"Longitude (Degrees)\")\n",
    "axs[2].set_ylabel(\"Latitude (Degrees)\")\n",
    "aspect_fig = axs[2].imshow(\n",
    "    read_aspect,\n",
    "    cmap=cmap_instance,\n",
    "    vmin=0,\n",
    "    origin=\"lower\",\n",
    "    vmax=360,\n",
    "    extent=(lat_min, lat_max, lon_min, lon_max),\n",
    ")\n",
    "\n",
    "\n",
    "cbar = fig.colorbar(\n",
    "    aspect_fig, fraction=0.046 * read_aspect.shape[0] / read_aspect.shape[1], pad=0.04\n",
    ")\n",
    "cbar_ticks = np.linspace(0, 360, 8)\n",
    "cbar.set_ticks(cbar_ticks)\n",
    "cbar.set_label(\"Aspect (Degrees)\")\n",
    "\n",
    "axs[3].set_xlim(lat_min, lat_max)\n",
    "axs[3].set_ylim(lon_min, lon_max)\n",
    "axs[3].set_title(\"Slope - TN 30m Resolution\")\n",
    "axs[3].set_xlabel(\"Longitude (Degrees)\")\n",
    "axs[3].set_ylabel(\"Latitude (Degrees)\")\n",
    "slope_fig = axs[3].imshow(\n",
    "    read_slope,\n",
    "    cmap=cmap_instance,\n",
    "    vmin=0,\n",
    "    origin=\"lower\",\n",
    "    vmax=65.9,\n",
    "    extent=(lat_min, lat_max, lon_min, lon_max),\n",
    ")\n",
    "\n",
    "\n",
    "cbar = fig.colorbar(\n",
    "    slope_fig, fraction=0.046 * read_slope.shape[0] / read_slope.shape[1], pad=0.04\n",
    ")\n",
    "cbar_ticks = np.linspace(0, 65.9, 8)\n",
    "cbar.set_ticks(cbar_ticks)\n",
    "cbar.set_label(\"Slope (Degrees)\")\n",
    "\n",
    "plt.subplots_adjust(wspace=0.4, hspace=0.6)\n",
    "plt.tight_layout()\n",
    "\n",
    "\n",
    "plt.show()\n",
    "\n",
    "print(\"You have successfully visualized the IDX files.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e80b3656-dc02-4849-9354-a9abf3f50c06",
   "metadata": {},
   "outputs": [],
   "source": [
    "PORT = \"8989\"  # Dont change this since this is the forwarded port.\n",
    "ADDRESS = \"0.0.0.0\"  # Local to the server; 0.0.0.0 or 127.0.0.1\n",
    "\n",
    "print(\"You have successfully set the port for the dashboard.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f57edc5-a3c5-420f-8989-a07996ca3148",
   "metadata": {},
   "outputs": [],
   "source": [
    "URL = \"https://maritime.sealstorage.io/api/v0/s3/utah/nsdf/somospie/tennessee_30m_terrain/Tennessee_terrain_parameters.idx?access_key=any&secret_key=any&endpoint_url=https://maritime.sealstorage.io/api/v0/s3&cached=arco\"\n",
    "print(\"You have successfully loaded the IDX Tennessee data from remote storage.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0596942a-9307-439d-bf27-e75c77a91ea0",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "## Stop this cell before running another visualization dashboard\n",
    "!python -m panel serve openvisuspy/src/openvisuspy/dashboards --log-file \"files/log.log\" --dev --allow-websocket-origin='*' --address=\"{ADDRESS}\" --port \"{PORT}\" --args \"{URL}\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f44e7138-bbbe-4368-8b2f-d6989830f7c0",
   "metadata": {},
   "source": [
    "##### Visit http://localhost:8989 to explore the dashboard"
   ]
  }
 ],
 "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.8.10"
  },
  "vscode": {
   "interpreter": {
    "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}