TerrainDashboardv2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterio\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "import hvplot.xarray\n",
    "import holoviews as hv\n",
    "import geoviews as gv\n",
    "import rioxarray as rx\n",
    "import os\n",
    "import panel as pn\n",
    "from holoviews import opts\n",
    "hv.extension('bokeh', width=90)\n",
    "pn.extension()\n",
    "import xrspatial as xrs\n",
    "from holoviews.streams import SingleTap\n",
    "from holoviews.streams import Tap\n",
    "from decimal import Decimal\n",
    "from math import isnan\n",
    "from math import ceil,floor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Creating variables which will store shared value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = xr.open_rasterio(\"https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/1/TIFF/n19w065/USGS_1_n19w065.tif\")[0,:,:].astype('float64')\n",
    "data.attrs = xr.open_rasterio(\"https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/1/TIFF/n19w065/USGS_1_n19w065.tif\").attrs\n",
    "data = data.where(data>0)\n",
    "maps = ['inferno','blues','kb','kgy','viridis']\n",
    "sloped = 0\n",
    "aspects = 0\n",
    "curvature = 0\n",
    "maxheight = 0\n",
    "minheight = 0\n",
    "maxslope = 0\n",
    "minslope = 0\n",
    "maxaspect = 0\n",
    "minaspect = 0\n",
    "mincurve = 0\n",
    "maxcurve = 0\n",
    "test = ''\n",
    "link = ''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Making useful helper functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def createlink(x,y):\n",
    "    \"\"\"\n",
    "    Creates the links which will be given to rasterio to pull the data. Returns a string\n",
    "    \"\"\"\n",
    "    global link\n",
    "    link = 'https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/1/TIFF/'\n",
    "    north = \"{:02d}\".format(round(abs(y)))\n",
    "    west =  \"{:03d}\".format(round(abs(x)))\n",
    "    coordinates = 'n'+north+'w'+west\n",
    "    finallink = link+coordinates+'/USGS_1_'+coordinates+'.tif'\n",
    "    link = finallink\n",
    "    return (finallink,coordinates) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getmaxmin(dset):\n",
    "    \"\"\"\n",
    "    calculates the maximum and minimum of an xarray data set and returns both in a tuple\n",
    "    \"\"\"\n",
    "    return (round(float(dset.min())),round(float(dset.max())))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getpercentage(orig,new):\n",
    "    \"\"\"\n",
    "    If given two xarray data sets (the original and the some augmented version of the orginal\n",
    "    which omits certian values) it will return what percent new is of the original\n",
    "    \"\"\"\n",
    "    origcount = int((~np.isnan(orig.values)).sum())\n",
    "    newcount = int((~np.isnan(new.values)).sum())\n",
    "    return round((newcount/origcount)*100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getregions(startpoint,endpoint):\n",
    "    \"\"\"\n",
    "    When given two bounding points the function will provide what regions are needed in order\n",
    "    represent all points in that region. The function will return a list containing tuples which\n",
    "    hold the the latititude and longitude of the points with the largest absoulte value in each\n",
    "    1 degree by 1 degree square.\n",
    "    \"\"\"\n",
    "    #NOTE: this only works if the startpoint is below the end point. Remember to compensate for that in code.\n",
    "    start = [floor(float(abs(startpoint[0]))),floor(float(abs(startpoint[1])))]\n",
    "    end = [ceil(float(abs(endpoint[0]))),ceil(float(abs(endpoint[1])))]\n",
    "    regions = []\n",
    "    for horiz in range(start[0]+1,end[0]+1):\n",
    "        for vert in range(start[1]+1,end[1]+1):\n",
    "            regions.append((horiz,vert))\n",
    "    return regions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getregiondata(points):\n",
    "    \"\"\"\n",
    "    When given a list from the getregions function, returns a list containing the xarray raster data for each\n",
    "    region\n",
    "    \"\"\"\n",
    "    regions = []\n",
    "    for p in points:\n",
    "        try:\n",
    "            place = createlink(p[0],p[1])[0]\n",
    "            data = rx.open_rasterio(place, masked=True).squeeze().astype('float64')\n",
    "            data.attrs = xr.open_rasterio(place).attrs\n",
    "            data.attrs['res'] = (1,np.cos(p[1]*(np.pi/180)))\n",
    "            regions.append(data)\n",
    "        except:\n",
    "            print('point '+str(p)+' is not avaliable')\n",
    "    return regions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bounddata(startpoint,endpoint):\n",
    "    \"\"\"\n",
    "    When given a start point and endpoint the function will return an aggregated raster with all the required 1 degree by 1 degree squares merged.\n",
    "    function must be given the corner of the region whose absolute values for the latitude and longitude are largest as the end point. Use the ensurepointscorrect function\n",
    "    to get around this. \n",
    "    \"\"\"\n",
    "    regions = getregions(startpoint,endpoint)\n",
    "    unprocessed = getregiondata(regions)\n",
    "    large = combineregions(unprocessed)\n",
    "    processed = large.sel(x=slice(endpoint[0],startpoint[0]),y=slice(endpoint[1],startpoint[1])).astype('float64')\n",
    "    processed.attrs['res'] = (1,1)\n",
    "    return processed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def combineregions(data):\n",
    "    \"\"\"\n",
    "    when givena list of xarray raster data returns merged version of all the data\n",
    "    \"\"\"\n",
    "    processed=[]\n",
    "    for tile in data:\n",
    "        # Removes 'halo region' which extands beyond the degree square\n",
    "        startpointx = round(float(tile.x[0]))\n",
    "        startpointy = round(float(tile.y[0]))\n",
    "        endpointx = round(float(tile.x[-1]))\n",
    "        endpointy = round(float(tile.y[-1]))\n",
    "        tile.name = 'box'\n",
    "        tile = tile.sel(x=slice(startpointx,endpointx),y=slice(startpointy,endpointy))\n",
    "        processed.append(tile)\n",
    "    return xr.combine_by_coords([sq.to_dataset() for sq in processed]).to_array()[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ensurepointscorrect(startpoint,endpoint):\n",
    "    \"\"\"\n",
    "    When given any two bounding points in tuples in latitude and longitude the function will return the corner\n",
    "    with the largest absolute value for the latitude and longitude as the endpoint and the point diaganally \n",
    "    opposite of it as the startpoint. Returns a tuple in the form (startpoint,endpoint) where start and end point\n",
    "    are also tuples\n",
    "    \"\"\"\n",
    "    vertcheck = startpoint[1]<endpoint[1]\n",
    "    horizcheck = startpoint[0]>endpoint[0]\n",
    "    if vertcheck and horizcheck:\n",
    "        return (startpoint,endpoint)\n",
    "    else:\n",
    "        dx = abs(startpoint[0]-endpoint[0])\n",
    "        dy = abs(startpoint[1]-endpoint[1])\n",
    "        starty = startpoint[1] if vertcheck else endpoint[1]\n",
    "        startx = startpoint[0] if horizcheck else endpoint[0]\n",
    "        endy = starty+dy\n",
    "        endx = startx-dx\n",
    "        return ((startx,starty),(endx,endy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Making widgets and Text Displays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmapwidget = pn.widgets.Select(name='Select Color Map', options=maps)\n",
    "heightquantileslider = pn.widgets.IntSlider(name='Height Quantile Slider', start=0, end=32, value=7)\n",
    "heightrangeslider = pn.widgets.RangeSlider(step=1,start = minheight, end = maxheight, value=(minheight,maxheight))\n",
    "sloperangeslider = pn.widgets.RangeSlider(start = minslope, end=maxslope,step=1,value=(minslope,maxslope))\n",
    "aspectrangeslider = pn.widgets.RangeSlider(start = minaspect, end=maxaspect,step=1,value=(minaspect,maxaspect))\n",
    "curvaturerangeslider = pn.widgets.RangeSlider(start = mincurve, end=maxcurve,step=1,value=(mincurve,maxcurve))\n",
    "selectlatp1 = pn.widgets.TextInput(name='Latitude1', value='19')\n",
    "selectlonp1 = pn.widgets.TextInput(name='Longitude1', value='-65')\n",
    "selectlatp2 = pn.widgets.TextInput(name='Latitude2', value='18')\n",
    "selectlonp2 = pn.widgets.TextInput(name='Longitude2', value='64')\n",
    "button = pn.widgets.Button(name='Process Points', button_type='primary')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "numclicksbefore = button.clicks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "heightp = pn.pane.Markdown()\n",
    "heightpercentpanel = pn.Column(heightp,css_classes=['panel-widget-box'],margin=25)\n",
    "slopep = pn.pane.Markdown()\n",
    "slopepercentpanel = pn.Column(slopep,css_classes=['panel-widget-box'],margin=25)\n",
    "aspectp = pn.pane.Markdown()\n",
    "aspectn = pn.pane.Markdown('# Aspect is the direction in which the land is pointing. The slider can be from 0-360 where 0 is absolute North')\n",
    "aspectpercentpanel = pn.Column(aspectp,css_classes=['panel-widget-box'],margin=25)\n",
    "aspectnote = pn.Column(aspectn,css_classes=['panel-widget-box'],margin=25)\n",
    "curvep = pn.pane.Markdown()\n",
    "curvaturepercentpanel = pn.Column(curvep,css_classes=['panel-widget-box'],margin=25)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Creating the displays for each tab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tiles1 = gv.tile_sources.EsriImagery().opts(width=600, height=550).redim.range(Latitude=(-50,50), Longitude=(-50,50))\n",
    "#tiles = gv.tile_sources.EsriImagery().opts(width=600, height=550).opts(xlim=(-50, 50), ylim=(-50, 50))\n",
    "\n",
    "tap_stream1 = SingleTap(source=tiles1, x=10, y=0)\n",
    "@pn.depends(tap_stream1.param.x, tap_stream1.param.y)\n",
    "def lat_lon1(x,y):\n",
    "    return f\"### Lat: {y:.4f} Lon: {x:.4f}\"\n",
    "\n",
    "def interactive_click1(x, y):\n",
    "    return gv.Points([(x,y)])\n",
    "\n",
    "click1 = hv.DynamicMap(interactive_click1, streams=[tap_stream1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#tiles2 = hv.element.tiles.EsriImagery().opts(width=600, height=550)\n",
    "\n",
    "tiles2=gv.tile_sources.EsriImagery().opts(width=600, height=550).redim.range(Latitude=(-50,50), Longitude=(-50,50))\n",
    "'''\n",
    "tap_stream2 = Tap(source=tiles2, x=10, y=0)\n",
    "\n",
    "@pn.depends(tap_stream2.param.x, tap_stream2.param.y)\n",
    "def lat_lon2(x,y):\n",
    "    return f\"### Lat: {y:.4f} Lon: {x:.4f}\"\n",
    "\n",
    "def interactive_click2(x, y):\n",
    "    return gv.Points([(x,y)])\n",
    "\n",
    "click2 = hv.DynamicMap(interactive_click2, streams=[tap_stream2])\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@pn.depends(colormap = cmapwidget, x=tap_stream1.param.x, y=tap_stream1.param.y)\n",
    "def clicklocation(colormap,x,y):\n",
    "    global data\n",
    "    global test\n",
    "    global sloped\n",
    "    global aspects\n",
    "    global curvature\n",
    "    try:\n",
    "        place = createlink(x,y)[0]\n",
    "        data = rx.open_rasterio(place, masked=True).squeeze().astype('float64')\n",
    "        data.attrs = xr.open_rasterio(place).attrs\n",
    "        data.attrs['res'] = (1,1)\n",
    "        data = data.where(data>0)\n",
    "        test = 'All good!'\n",
    "    except:\n",
    "        test = 'There was an error'\n",
    "    sloped = xrs.slope(data)\n",
    "    sloped.attrs['res']=(1,1)\n",
    "    aspects = xrs.aspect(data)\n",
    "    aspects.attrs['res']=(1,1)\n",
    "    curvature = xrs.slope(sloped)\n",
    "    curvature.attrs['res']=(1,1)\n",
    "    (heightrangeslider.start,heightrangeslider.end)=getmaxmin(data)\n",
    "    heightrangeslider.value = getmaxmin(data)\n",
    "    (sloperangeslider.start,sloperangeslider.end)=getmaxmin(sloped)\n",
    "    sloperangeslider.value = getmaxmin(sloped)\n",
    "    (aspectrangeslider.start,aspectrangeslider.end)=getmaxmin(aspects)\n",
    "    aspectrangeslider.value = getmaxmin(aspects)\n",
    "    (curvaturerangeslider.start,curvaturerangeslider.end)=getmaxmin(curvature)\n",
    "    curvaturerangeslider.value = getmaxmin(curvature)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@pn.depends(colormap = cmapwidget, clicker = button.param.clicks)\n",
    "def inputlocation(colormap,clicker):\n",
    "    global data\n",
    "    global test\n",
    "    global sloped\n",
    "    global aspects\n",
    "    global curvature\n",
    "    global numclicksbefore\n",
    "    global sample\n",
    "    numclicksbefore = clicker\n",
    "    points = ensurepointscorrect((float(selectlonp1.value),float(selectlatp1.value)),(float(selectlonp2.value),float(selectlatp2.value)))\n",
    "    sample = selectlonp1.value\n",
    "    data = bounddata(points[0],points[1])\n",
    "    data = data.where(data>0)\n",
    "    test = 'All good!'\n",
    "    #except Exception as e: test = e\n",
    "    sloped = xrs.slope(data)\n",
    "    aspects = xrs.aspect(data)\n",
    "    curvature = xrs.curvature(data)\n",
    "    (heightrangeslider.start,heightrangeslider.end)=getmaxmin(data)\n",
    "    heightrangeslider.value = getmaxmin(data)\n",
    "    (sloperangeslider.start,sloperangeslider.end)=getmaxmin(data)\n",
    "    sloperangeslider.value = getmaxmin(data)\n",
    "    (aspectrangeslider.start,aspectrangeslider.end)=getmaxmin(data)\n",
    "    aspectrangeslider.value = getmaxmin(data)\n",
    "    (curvaturerangeslider.start,curvaturerangeslider.end)=getmaxmin(data)\n",
    "    curvaturerangeslider.value = getmaxmin(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@pn.depends(x=tap_stream1.param.x, y=tap_stream1.param.y, clicker = button.param.clicks)\n",
    "def shaded(x,y,clicker):\n",
    "    global data\n",
    "    shade = xrs.hillshade(data)\n",
    "    return shade.hvplot.image(x = 'x', y = 'y', rasterize = True, geo = True, tiles = 'OSM', cmap = cmapwidget, alpha = 20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@pn.depends(x=tap_stream1.param.x, y=tap_stream1.param.y, level = heightquantileslider, clicker = button.param.clicks)\n",
    "def quantile(x,y, level,clicker):\n",
    "    global data\n",
    "    quant = xrs.quantile(data,k=level,ignore_vals=(0,))\n",
    "    return quant.hvplot.image(x = 'x', y = 'y',rasterize = True, geo = True, tiles = 'OSM', cmap = cmapwidget, alpha = 20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@pn.depends(x=tap_stream1.param.x, y=tap_stream1.param.y, border = heightrangeslider, clicker = button.param.clicks)\n",
    "def heightrange(x,y, border, clicker):\n",
    "    global data\n",
    "    manipulate = data.where(np.logical_and(data>border[0],data<border[1]))\n",
    "    heightp.object = '# Your selected range is {0} % of this region'.format(str(getpercentage(data,manipulate)))\n",
    "    shade = xrs.hillshade(manipulate)    \n",
    "    return shade.hvplot.image(x = 'x', y = 'y',rasterize = True, geo = True, tiles = 'OSM', cmap = cmapwidget)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@pn.depends(x=tap_stream1.param.x, y=tap_stream1.param.y, border = sloperangeslider, clicker = button.param.clicks)\n",
    "def slope(x,y, border,clicker):\n",
    "    global sloped\n",
    "    global tester\n",
    "    manipulate = sloped.where(np.logical_and(sloped>border[0],sloped<border[1]))\n",
    "    tester=border\n",
    "    slopep.object = '# Your selected range is '+ str(getpercentage(sloped,manipulate))+'% of this region'\n",
    "    #slopep.object = '# Your selected range is {0} % of this region'.format(str(getpercentage(sloped,manipulate)))\n",
    "    return manipulate.hvplot.image(x = 'x', y = 'y',rasterize = True, geo = True, tiles = 'OSM', cmap = cmapwidget)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@pn.depends(x=tap_stream1.param.x, y=tap_stream1.param.y, border = aspectrangeslider, clicker = button.param.clicks)\n",
    "def aspect(x,y, border,clicker):\n",
    "    global aspect\n",
    "    manipulate = aspects.where(np.logical_and(aspects>border[0],aspects<border[1]))\n",
    "    aspectp.object = '# Your selected range is '+ str(getpercentage(aspects,manipulate))+'% of this region'\n",
    "    #aspectp.object = '# Your selected range is {0} % of this region'.format(str(getpercentage(aspect,manipulate)))\n",
    "    return manipulate.hvplot.image(x = 'x', y = 'y',rasterize = True, geo = True, tiles = 'OSM', cmap = cmapwidget)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@pn.depends(x=tap_stream1.param.x, y=tap_stream1.param.y, border = curvaturerangeslider,clicker = button.param.clicks)\n",
    "def curvature(x,y, border,clicker):\n",
    "    global curvature\n",
    "    manipulate = curvature.where(np.logical_and(curvature>border[0],curvature<border[1]))\n",
    "    curvep.object = '# Your selected range is '+ str(getpercentage(curvature,manipulate))+'% of this region'\n",
    "    #curvep.object = '# Your selected range is {0} % of this region'.format(str(getpercentage(curvature,manipulate)))\n",
    "    return curvature.hvplot.image(x = 'x', y = 'y',rasterize = True, geo = True, tiles = 'OSM', cmap = cmapwidget)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Creating the canvas for each tab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "selector1 = pn.Column(pn.panel(lat_lon1, width=500), (tiles1 * click1))\n",
    "#selector2 = pn.Column(pn.panel(lat_lon2, width=500), (tiles2 * click2))\n",
    "clicker = pn.Column(cmapwidget,clicklocation,selector1)\n",
    "points = pn.Column(pn.Row(selectlatp1,selectlonp1),\n",
    "                   pn.Row(selectlatp2,selectlonp2),\n",
    "                   button\n",
    "                  )\n",
    "mappart = pn.Column(tiles2,cmapwidget)\n",
    "enterlocation = pn.Row(mappart,points,inputlocation)\n",
    "quantileview = pn.Row(quantile, heightquantileslider)\n",
    "heightrangeview = pn.Row(heightrange, pn.Column(heightrangeslider,heightpercentpanel))\n",
    "slopeview = pn.Row(slope, pn.Column(sloperangeslider,slopepercentpanel))\n",
    "aspectview = pn.Row(aspect, pn.Column(aspectrangeslider,aspectpercentpanel),aspectnote)\n",
    "curvatureview = pn.Row(curvature, pn.Column(curvaturerangeslider,curvaturepercentpanel))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "option = pn.Tabs(('Location Clicker',clicker),\n",
    "                 ('Input Location',enterlocation)\n",
    "                )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#thing = pn.Row(shaded, button.on_click(shaded))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dashboard = pn.Tabs(('Options',option), \n",
    "                    (\"Shaded View\", shaded),\n",
    "                    ('Height Range', heightrangeview),\n",
    "                    ('Height Quantile', quantileview),\n",
    "                    ('Aspect View', aspectview),\n",
    "                    ('Slope View', slopeview),\n",
    "                    ('Curvature View', curvatureview),\n",
    "                    dynamic = True\n",
    "                   )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dashboard.servable()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:notebook] *",
   "language": "python",
   "name": "conda-env-notebook-py"
  },
  "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.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}