TouchTerrain_standalone_jupyter_notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TouchTerrain standalone in a jupyter notebook\n",
    "Chris Harding, Dec. 24, 2022 (<charding@iastate.edu>)\n",
    "\n",
    "- This jupyter notebook runs a standalone version of TouchTerrain, similar to `TouchTerrain_standalone.py`\n",
    "- This notebook needs to be run in Python __3.x__  \n",
    "- If you're using this via the touchterrain_jupyter docker container, ignore all the install stuff that follows. The image will already contain nearly all packages, except for the touchterrain module itself, which you should already have installed via the install_touchterrain.sh shell script. You might want to rename/save this notebook into something other than TouchTerrain_standalone_jupyter_notebook.ipnb b/c if you ever update the touchterrain module from github (update_touchterrain.sh) it will overwrite it!\n",
    "\n",
    "\n",
    "- Use pip to build and install a package called touchterrain (all lowercase!), which contains all functions needed to run in standalone and in server mode. This will also install all dependencies needed(!)\n",
    "- you can run pip either in a terminal or inside a jupyter cell (see cell below)\n",
    "- to run it in a separate terminal, `cd` to the folder that contains the `setup.py` file (which will contain a folder called touchterrain)\n",
    "    - enter:  `pip install .` \n",
    "    - note the . at the end, which will make pip run the script in setup.py\n",
    "    - this will install all dependencies (including earthengine-api) \n",
    "- You will also need to install __gdal__ which can be tricky when using pip as it requires a C++ compiler. Because of many Windows users not wanting to go that route, it's commented out in setup.py, but if you do have a working compiler you can uncomment it and see if it works with pip.\n",
    "- On Windows you could get the whl file for your version of Python from [here](https://www.lfd.uci.edu/~gohlke/pythonlibs/) and install it separately with `pip install <whl file>` (Note that this is sadly no longer being updated!)\n",
    "- gdal is also available via [conda-forge](https://anaconda.org/conda-forge/gdal) and if you have conda installed you could try: `conda install -c conda-forge gdal` \n",
    "- Finally. there are two optional modules which you may want to install via pip as well:\n",
    "    - `pip install geemap`  will allow you to use an interactive in-cell map to digitize the print area\n",
    "    - `pip install k3d` will allow you to preview your STL file in a 3D viewer\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# if you have already installed touchterrain via your own terminal, or if you are using the docker version, skip this cell!\n",
    "\n",
    "# Uncomment these 2 lines and run this cell (Shift-Enter) to have pip install the touchterrain module and \n",
    "# all its dependencies. Warning: this can take a a while! You may not see any message during the install \n",
    "# but you should see the full log when pip is done (or when you got errors)\n",
    "# You'll need something similar to install then optional modules mentioned above.\n",
    "\n",
    "#import sys\n",
    "#!{sys.executable} -m pip install . "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### If you're not going to use Earth Engine's (EE) online DEM data, skip the next section and go to _Running Touchterrain standalone_"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preparing to use Google Earth Engine online DEM rasters \n",
    "\n",
    "- You don't need this if you only want to import terrain from locally stored raster files (geotiffs)!\n",
    "- TouchTerrain can use DEM data from Google Earth Engine (given the corners of the area), but you need to first go to https://signup.earthengine.google.com/ This should work with a standard Google account, with which you later need to authenticate with Earth Engine to get access to online DEM data\n",
    "\n",
    "- The authentication requires that you have this authentication file in place: `<your home folder/.config/earthengine/credentials`\n",
    "- This credentials file is created by running `ee.Authenticate()` __once__(!). \n",
    "- To create the credentials file uncomment `ee.Authenticate()` in the next cell and run it once. \n",
    "- After you have done this ONCE, you should not have to use `ee.Authenticate()` anymore, so comment it out again. (`#ee.Authenticate()`)\n",
    "- You will still need to run this cell for its `import ee` every time, just not the `ee.Authenticate()` part\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ee # earthengine-api should have been installed via pip earlier\n",
    "#ee.Authenticate()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- When you run: `ee.Authenticate()` this text should appear and you should be redirected to a webpage:\n",
    "- (Note: you may want to open this link in a new browser tab)\n",
    "```\n",
    "To authorize access needed by Earth Engine, open the following URL in a web browser and follow the instructions. If the web browser does  not start automatically, please manually browse the URL below.\n",
    "\n",
    "    https://accounts.google.com/o/oauth2/auth?client_id=5172225062...\n",
    "\n",
    "```\n",
    "- the web page will have you select a Google account for use with ee and give you a code.\n",
    "- For more info, see `EarthEngine_authentication_guide.md`\n",
    "- paste in your code here: `Enter verification code: <your code>` and hit Enter - you should get `Successfully saved authorization token.`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running Touchterrain standalone (local and online DEM data)\n",
    "\n",
    "- Put your settings into the dictionary below and hit __Shift-Enter__\n",
    "- For more info on the settings, look at the ReadMe on https://github.com/ChHarding/TouchTerrain_for_CAGEO\n",
    "- Note, however, that the settings given below are in Python syntax, whereas the ReadMe describes the JSON syntax used in the config file\n",
    "- Both are very similar, except for None and True/False\n",
    "\n",
    "``` \n",
    "    Python:  JSON:\n",
    "    None     null\n",
    "    True     true\n",
    "    False    false\n",
    "```\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# import packages\n",
    "import os, sys\n",
    "from pprint import pprint\n",
    "\n",
    "# The touchterrain module should have been installed via pip earlier\n",
    "# This will also run ee.Initialize() and should show: EE init() worked with .config/earthengine/credentials\n",
    "# If you don't plan on using online DEM data and thus didn't authenticate earlier, you will get a warning, which you can ignore.\n",
    "from touchterrain.common import TouchTerrainEarthEngine as TouchTerrain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "args = {\n",
    "    # DEM/Area to print\n",
    "    \n",
    "    # A: use local DEM raster (geotiff)\n",
    "    #\"importedDEM\": \"stuff/pyramid.tif\",  # path to the geotif in relation to where this notebook sits\n",
    "    \n",
    "    # B: use area and a DEM online source via EarthEngine\n",
    "    \"importedDEM\": None,\n",
    "    \"DEM_name\": \"USGS/3DEP/10m\",   # DEM source\n",
    "    # the following defines the area, but you can also define it by hand (see digitizing) \n",
    "    \"bllat\": 44.50185267072875,   # bottom left corner lat\n",
    "    \"bllon\": -108.25427910156247, # bottom left corner long\n",
    "    \"trlat\": 44.69741706507476,   # top right corner lat\n",
    "    \"trlon\": -107.97962089843747, # top right corner long\n",
    "\n",
    "    # 3D print parameters\n",
    "    \"tilewidth\": 120,  # width of each tile in mm, (tile height will be auto calculated)\n",
    "    \"printres\": 0.4,  # resolution (horizontal) of 3D printer in mm, should be your NOZZLE size or just a bit less! \n",
    "                      # Using something like 0.01 will NOT print out a super detailed version as you slicer will remove\n",
    "                      # super fine details anyway! You'll just wait a long time and get a super large STL file!\n",
    "    \n",
    "    \"ntilesx\": 1, # number of tiles in x  \n",
    "    \"ntilesy\": 1, # number of tiles in y    \n",
    "\n",
    "    \"basethick\": 0.6,   # thickness (in mm) of printed base\n",
    "    \"zscale\": 2,      # elevation (vertical) scaling\n",
    "    \"fileformat\": \"STLb\",  # format of 3D model files: \"obj\" wavefront obj (ascii),\n",
    "                           #   \"STLa\" ascii STL or \"STLb\" binary STL.\n",
    "                           #   To export just the (untiled) raster (no mesh), use \"GeoTiff\" \n",
    "    \"zip_file_name\": \"myterrain\",   # base name of zipfile, .zip will be added\n",
    "\n",
    "    # Expert settings\n",
    "    \"tile_centered\": False, # True-> all tiles are centered around 0/0, False, all tiles \"fit together\"\n",
    "    \"CPU_cores_to_use\" : 0, # 0: use all available cores, None: don't use multiprocessing (single core only)\n",
    "                            # multi-core will be much faster for more than 1 tile \n",
    "    \"max_cells_for_memory_only\" : 5000^2, # if number of raster cells is bigger than this, use temp_files instead of memory.\n",
    "                            # set this very high to force use of memory and lower it if you run out of memory\n",
    "    \"no_bottom\": False,   # omit bottom triangles? Most slicers still work and it makes smaller files\n",
    "    \"no_normal\": True,    # Don't calculate normals for triangles. This is significantly faster but some 3D model viewers may need them.\n",
    "    \"bottom_image\": None, # \"stuff/TouchTerrain_bottom_example.png\", # 1 band greyscale image used for bottom relief\n",
    "    \"ignore_leq\": None,   # set all values <= this to NaN so they don't print\n",
    "    \"lower_leq\": None,    # e.g. [0.0, 2.0] values <= 0.0 will be lowered by 2mm in the final model\n",
    "    \"unprojected\": False, # don't project to UTM (for EE rasters only)\n",
    "    \"projection\": None,   # None means use the closest UTM zone. Can be a EPSG number (int!) instead but not all work. \n",
    "    \"only\" : None,        # if not None: list with x and y tile index (1 based) of the only tile to process\n",
    "                          #   e.g. [1,1] will only process the tile in upper left corner, [2,1] the tile right to it, etc.\n",
    "    \"importedGPX\": [],    # list of gpx path file(s) to be use (optional: see next cell)\n",
    "}\n",
    "########################################################\n",
    "\n",
    "# if we want to work on a local raster, get the full pathname to it\n",
    "if args[\"importedDEM\"] != None: \n",
    "    from os.path import abspath\n",
    "    args[\"importedDEM\"]= abspath(args[\"importedDEM\"]) \n",
    "    print(\"reading in local DEM:\", args[\"importedDEM\"])\n",
    "print(\"settings stored, ready to process\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optional: Drape GPX path file(s) over the terrain (thanks to KohlhardtC!)\n",
    "- you can drape one or more gpx (path) files over your terrain\n",
    "- the gpx file(s) have to be part of list (`importedGPX`)\n",
    "- `gpxPathHeight` (meters) defines how much a path is elevated above the terrain. Use a negative number to create a trench.\n",
    "- `gpxPixelsBetweenPoints` (meters) lets you reduce the number of point in your path to place a point only every X meters. This can help to simplify complex patterns.\n",
    "- `gpxPathThickness` (meters) controls the thickness of the path\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Note: you must uncomment the last line in this cell to actually use these gpx settings!\n",
    "from touchterrain.common.TouchTerrainGPX import *\n",
    "gpx_args = {   \n",
    "    # Area for using the example GPX test files\n",
    "    \"bllat\": 39.32205105794382,   # bottom left corner lat\n",
    "    \"bllon\": -120.37497608519418, # bottom left corner long\n",
    "    \"trlat\": 39.45763749030933,   # top right corner lat\n",
    "    \"trlon\": -120.2002248034559, # top right corner long\n",
    "\n",
    "    \"importedGPX\": # gpx example files.\n",
    "                 [\"stuff/gpx-test/DLRTnML.gpx\",\n",
    "                  \"stuff/gpx-test/DonnerToFrog.gpx\",\n",
    "                  \"stuff/gpx-test/CinTwistToFrog.gpx\",\n",
    "                  \"stuff/gpx-test/sagehen.gpx\",\n",
    "                  \"stuff/gpx-test/dd-to-prosser.gpx\",\n",
    "                  \"stuff/gpx-test/alder-creek-to-crabtree-canyon.gpx\",\n",
    "                  \"stuff/gpx-test/ugly-pop-without-solvang.gpx\",  \n",
    "                 ],\n",
    "     \"gpxPathHeight\": 10,  # Currently we plot the GPX path by simply adjusting the\n",
    "                           # raster elevation at the specified lat/lon,\n",
    "                           # therefore this is in meters. Negative numbers are ok \n",
    "                           # and put a dent in the mdoel  \n",
    "     \"gpxPixelsBetweenPoints\" : 20, # GPX Files haves a lot of points. A higher \n",
    "                                    #number will create more space between lines drawn\n",
    "                                    # on the model and can have the effect of making the paths look a bit cleaner \n",
    "     \"gpxPathThickness\" : 2, # Stack paralell lines on either side of primary line \n",
    "                             # to create thickness. A setting of 1 probably looks the best\n",
    "}\n",
    "\n",
    "# uncomment the next line if you want to use gpx_args!\n",
    "#args = {**args, **gpx_args}; print(args) # merge gpx_args into args, ** unrolls dicts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optional: Show map and digitize print area shape (box, circle or polygon)\n",
    "- the next 2 cells will show you a hillshaded Google map and will let you manually define a shape, either a rectangle (box), a circle or a polygon\n",
    "- it will show you a red box of the area defined earlier (with bllat, etc.) in case you want to digitize inside it. However, you're free to digitize anywhere on the map! Your digitized shape will always override the red box.\n",
    "- This requires the `geemap` module. To install it, type `pip install geemap` into a OS terminal (Again, not needed if you're using Docker)\n",
    "- Run the next cell to see the Google map\n",
    "- On the left, use the __Draw a Rectangle__, __Draw a Circle__ or __Draw a Polygon__ to define the shape of your print area\n",
    "- Run the cell below the Map to make a GeoJSON polygon and add it to args. (If you don't want to use the shape (or didn't digitize any) and want to use the red bounding box, just skip this step.)\n",
    "- Run the processing cell "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Geemap\n",
    "import geemap\n",
    "\n",
    "# Create an interactive map and center on default area\n",
    "center_lat = (args[\"trlat\"] + args[\"bllat\"]) / 2\n",
    "center_lon = (args[\"trlon\"] + args[\"bllon\"]) / 2\n",
    "Map = geemap.Map(center=(center_lat, center_lon), zoom=7) \n",
    "\n",
    "# make a hillshade layer and add it to map\n",
    "dem = ee.Image(args[\"DEM_name\"]) # DEM source\n",
    "hs = ee.Terrain.hillshade(dem, 315, 35) # sun azimuth and angle \n",
    "vis_params = {'min': 0,'max': 255,} # greyscale color ramp\n",
    "Map.addLayer(hs, vis_params, 'hillshade', shown=True, opacity=0.5) # semi transparent overlay\n",
    "\n",
    "# if GPX files were used, add them to the map\n",
    "if args[\"importedGPX\"] != None and len(args[\"importedGPX\"]) > 0: \n",
    "        gpx = ee.Feature(convert_to_GeoJSON(args[\"importedGPX\"]))\n",
    "        Map.addLayer(gpx, {\"color\":'00FFFF', \"strokeWidth\":\"1\"}, \"GPX line\", opacity=0.9)\n",
    "\n",
    "# show the currently defined bounding box\n",
    "# I'm commenting this out b/c I think geemap users don't really need to see that\n",
    "# given that they probably use geemap to interactively select that box. But, if you do\n",
    "# want to see it for some reason, uncomment the block below \n",
    "#rect = ee.Geometry.Rectangle(\n",
    "#        args[\"bllon\"], args[\"bllat\"], args[\"trlon\"], args[\"trlat\"])\n",
    "#rect_feature = ee.Feature(rect)\n",
    "#Map.addLayer(rect_feature, {\"color\": 'FF0000'}, \"bounding box\", opacity=0.5)\n",
    "\n",
    "Map # makes the interactive map show up as output of this cell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# in the map GUI above, digitize a rectangle, circle or polygon, then run this cell\n",
    "# (if you screw up, just digitize a new shape and run this cell again, it will always use\n",
    "# the last digitized shape ...)\n",
    "\n",
    "# make GeoJSON polygon from (last) digitized polygon feature\n",
    "polyft = Map.draw_last_feature\n",
    "from geojson import Polygon\n",
    "coords = polyft.getInfo()['geometry']['coordinates']\n",
    "poly = Polygon(coords) # coords[0] is the digitized polygon, [1] etc. would be holes\n",
    "args[\"polygon\"] = poly # use this GeoJSON polygon in processing \n",
    "print(\"Using digitized polygon with\", len(coords[0]), \"points\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Processing\n",
    "- Running the cell below processes the data and creates a zip file with the model file(s) inside. This zip file will be inside the tmp folder (which is inside the same folder your notebook file is in).\n",
    "- This may take some time! \n",
    "- During processing you'll see a star indicator (`In[*]`) and some log messages. (Those messages will also be in the logfile inside the zip)\n",
    "- You may see some red messages with 10%, etc. - don't worry, that's normal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "\n",
    "totalsize, full_zip_file_name = TouchTerrain.get_zipped_tiles(**args) # args are in a dict\n",
    "print(\"\\nDONE!\\n\\nCreated zip file\", full_zip_file_name,  \"%.2f\" % totalsize, \"Mb\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "# If you want to unzip the zip file, run this cell\n",
    "# (You will need to do this before using k3d for visualization)\n",
    "\n",
    "import os.path\n",
    "from glob import glob\n",
    "folder, file = os.path.splitext(full_zip_file_name) # get folder of zip file\n",
    "\n",
    "# unzip the zipfile into the folder it's already in\n",
    "import zipfile\n",
    "zip_ref = zipfile.ZipFile(full_zip_file_name, 'r')\n",
    "zip_ref.extractall(folder)\n",
    "zip_ref.close()\n",
    "print (\"unzipped files from\", full_zip_file_name, \"into the folder\", folder)\n",
    "print (folder, \"contains these files:\")\n",
    "for f in glob(folder + os.sep + \"*.*\"): print(\" \", f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize the STL file(s)\n",
    "\n",
    "- If you want to visualize your model(s), install k3d (`pip install k3d`) and run the cell below\n",
    "- this won't work with OBJ files, as k3d can't read them in.\n",
    "- binary STL should work, but ascii may not.\n",
    "- each tile will have a random color, so you can see the individual tiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# install k3d using the jupyter magic %pip command\n",
    "%pip install k3d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If you're running this notebook in Docker, you need to also install\n",
    "# and activate the jupyter k3d extension\n",
    "!jupyter nbextension install --py --user k3d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# activate extension\n",
    "!jupyter nbextension enable k3d --user --py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import k3d\n",
    "from glob import glob\n",
    "\n",
    "# get all stl files (tiles) in that folder\n",
    "mesh_files = glob(folder + os.sep + \"*.STL\")\n",
    "\n",
    "plot = k3d.plot()\n",
    "\n",
    "# Add all tiles with a random color\n",
    "from random import randint\n",
    "for m in mesh_files:\n",
    "    col = (randint(0,255) << 16) + (randint(0,255) << 8) + randint(0,255) # random rgb color as hex\n",
    "    print(\"adding to viewer:\", m, hex(col))\n",
    "    buf = open(m, 'rb').read()\n",
    "    plot += k3d.stl(buf, color=col)\n",
    "plot.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If you plan to refine your parameters and create different models, you should run this cell to remove \n",
    "# the preview folder first, otherwise the preview could show you old models.\n",
    "\n",
    "from shutil import rmtree\n",
    "rmtree(folder)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.6 64-bit",
   "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.6"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  },
  "vscode": {
   "interpreter": {
    "hash": "b1be0222797457612977eb3eaa00946922d56583fa52c09202534133b3ce7849"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}