WMTS_Stitcher_(V4).ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from holoviews.streams import PointDraw, PolyEdit, BoxEdit, PolyDraw\n",
    "from requests.packages.urllib3.util.retry import Retry\n",
    "from IPython.display import display, clear_output\n",
    "from ipywidgets import widgets, HBox, Label\n",
    "from requests.adapters import HTTPAdapter\n",
    "from ipywidgets import widgets, Layout\n",
    "from rasterio.io import MemoryFile\n",
    "from itertools import product\n",
    "from io import BytesIO\n",
    "from PIL import Image\n",
    "import cartopy.crs as ccrs\n",
    "import holoviews as hv\n",
    "import geoviews as gv\n",
    "import numpy as np\n",
    "import requests\n",
    "import rasterio\n",
    "import warnings\n",
    "import math\n",
    "import sys\n",
    "import os\n",
    "\n",
    "hv.extension('bokeh')\n",
    "warnings.filterwarnings('ignore')\n",
    "# Helpful URLS:\n",
    "    # https://stackoverflow.com/questions/15431044/can-i-set-max-retries-for-requests-request"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "url_list = ['http://c.tile.openstreetmap.org/{Z}/{X}/{Y}.png',\n",
    "            'https://s.basemaps.cartocdn.com/light_all/{Z}/{X}/{Y}.png',\n",
    "            'http://tile.stamen.com/terrain/{Z}/{X}/{Y}.png',\n",
    "            'https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{Z}/{Y}/{X}',\n",
    "            'https://server.arcgisonline.com/ArcGIS/rest/services/World_Topo_Map/MapServer/tile/{Z}/{Y}/{X}',\n",
    "            'https://server.arcgisonline.com/ArcGIS/rest/services/World_Terrain_Base/MapServer/tile/{Z}/{Y}/{X}',\n",
    "            'https://server.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer/tile/{Z}/{Y}/{X}',\n",
    "            'https://server.arcgisonline.com/ArcGIS/rest/services/World_Shaded_Relief/MapServer/tile/{Z}/{Y}/{X}',\n",
    "            'https://server.arcgisonline.com/ArcGIS/rest/services/World_Physical_Map/MapServer/tile/{Z}/{Y}/{X}',\n",
    "            'https://server.arcgisonline.com/ArcGIS/rest/services/USA_Topo_Maps/MapServer/tile/{Z}/{Y}/{X}',\n",
    "            'https://server.arcgisonline.com/ArcGIS/rest/services/Ocean_Basemap/MapServer/tile/{Z}/{Y}/{X}',\n",
    "            'https://server.arcgisonline.com/ArcGIS/rest/services/NatGeo_World_Map/MapServer/tile/{Z}/{Y}/{X}']\n",
    "\n",
    "url_selection = widgets.Dropdown(options=url_list, description='Tile URL:', disabled=False, \n",
    "                              layout=Layout(width='75%', height='40px'))\n",
    "display(url_selection)\n",
    "\n",
    "zoom_widget_value = widgets.IntSlider(min=0,max=19,step=1,description='Zoom Level:',orientation='horizontal')\n",
    "display(zoom_widget_value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "url = url_selection.value\n",
    "tiles = gv.WMTS(url, extents=(-180, -90, 180, 90), crs=ccrs.PlateCarree()).options(width=500, height=500)\n",
    "boxes = gv.Polygons([], crs=ccrs.PlateCarree()).options(fill_alpha=0.4, color='blue', line_color='blue', line_width=2)\n",
    "box_stream = BoxEdit(source=boxes)\n",
    "tiles * boxes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%time\n",
    "xs, ys = box_stream.element.array().T\n",
    "bbox = list(gv.util.project_extents((xs[0], ys[0], xs[2], ys[1]), ccrs.GOOGLE_MERCATOR, ccrs.PlateCarree()))\n",
    "print(\"Boudning box: \", bbox)\n",
    "\n",
    "zoom_level = zoom_widget_value.value\n",
    "# Get the basemap tiles for the given coordinates.\n",
    "convereted_lat_tile_list = [int(1/(2*math.pi)*2**zoom_level*(math.pi-math.log(math.tan(math.pi\n",
    "                                                /4+math.radians(y)/2)))) for y in [bbox[3], bbox[1]]]\n",
    "convereted_long_tile_list = [int(2**(zoom_level-1)*(x/180+1)) for x in [bbox[0], bbox[2]]]\n",
    "convereted_lat_tile_list.sort()\n",
    "convereted_long_tile_list.sort()\n",
    "\n",
    "# Get the top left and bottom right hand corner coordinates of where the actual basemap tile starts.\n",
    "# The top left hand corner of the bbox will not have the coordinates as the basemap tile because of \n",
    "# - them never going to line up perfectly.\n",
    "# We have to increment the long and lat of the bottom right hand corner because we need to get the next\n",
    "# - tile and caculate the coordinates of its top left hand corner.\n",
    "# Why I am doing this is because we have to geo-reference the image and have the exact coordinates of the\n",
    "# - actual image that we will stitch together from these tiles.\n",
    "first_lat = math.degrees(2*(math.atan(math.exp(math.pi-(convereted_long_tile_list[0]*2*math.pi)\n",
    "                                               /2**zoom_level))-math.pi/4))\n",
    "second_lat = math.degrees(2*(math.atan(math.exp(math.pi-((convereted_long_tile_list[-1] + 1)*2*math.pi)\n",
    "                                                /2**zoom_level))-math.pi/4))\n",
    "\n",
    "first_long = math.degrees(convereted_lat_tile_list[0]*2*math.pi/(2**zoom_level) - math.pi)\n",
    "second_long = math.degrees((convereted_lat_tile_list[-1] + 1)*2*math.pi/(2**zoom_level) - math.pi)\n",
    "\n",
    "# Have to create a range item so that we can do a product on it later in the algorithm. So lets increment,\n",
    "# - the last element in either list to include up to the next number. If we didn't then it would stop short\n",
    "# - one number.\n",
    "# Next I have to figure out the number of tiles for x and y, so that I can see how big the image is going \n",
    "# - to be.\n",
    "# We have to increment the last elemeent in the list because Python range function is not inclusive of the last\n",
    "# - number.\n",
    "lat_range_list = range(convereted_lat_tile_list[0], (convereted_lat_tile_list[1] + 1))\n",
    "number_of_y_tiles = len(lat_range_list)\n",
    "\n",
    "long_range_list = range(convereted_long_tile_list[0], (convereted_long_tile_list[1] + 1))\n",
    "number_of_x_tiles = len(long_range_list)\n",
    "\n",
    "total_number_of_tiles = number_of_x_tiles * number_of_y_tiles\n",
    "\n",
    "if (number_of_x_tiles * number_of_y_tiles) > 1000:\n",
    "    print(\"This might take a long time...\")\n",
    "print(\"Total number of images to get: \", total_number_of_tiles)\n",
    "\n",
    "# This loop will create the image from the previous calculated size, by multiplying the number of tiles for\n",
    "# - x and y by 256 (pixels). This is the standard size for every tile (256,256) for open street map. Next, \n",
    "# - I will do the product on the range of both the lat and long. This will give me all the tiles from the \n",
    "# - top left corner to the the bottom right. Then, I have a dictionary that sets the z, x, and y to the upper\n",
    "# - and lower value. I do this because there will be an Z, X, and Y at the end of the URL that I will want to inject\n",
    "# - the values into, and won't know if the incoming x, y, z will be capitialized or lowercase. Next, I format \n",
    "# - the url with the given zoom_level, and calculated tile number to download. Next, I setup a request session to make\n",
    "# - sure that I can have retries to download a tile if one fails. Then, If I get a good response then I check \n",
    "# - to see if we are on the first iteration, to set my min x and y tile numbers. I do this because I need the \n",
    "# - baseline to subtract all future tile numbers from in order to acurately know where each tile will go in the \n",
    "# - grid. This has to be done this way because not all tiles will be retrieved in order, and this saves lines \n",
    "# - of code by not having to make sure that the pictures are all in order. Next, I increment my iterator so I don't \n",
    "# - come back into the if statment, I open the image in memory, then subtract the baseline x and y from the current \n",
    "# - tile number, and multiply by the number of pixels. Ex: (4 - 2) * 256 = 512 for x. This will put the picture right \n",
    "# - at the calculated number of pixels either over or down. Next I stitch the opened image to the one I created\n",
    "# - at the beinging of the loop. This will hold all the images once they are all stitched in. Lastly, I close\n",
    "# - the image that is opened in memory.\n",
    "loop_iterator = 0\n",
    "stitched_image = Image.new('RGB', (number_of_x_tiles * 256, number_of_y_tiles * 256 ))\n",
    "for x, y in product(long_range_list, lat_range_list):\n",
    "    key = {\"Z\": zoom_level, \"z\": zoom_level, \"X\": x, \"x\": x, \"Y\": y, \"y\": y}\n",
    "    \n",
    "    internet_status = requests.Session()\n",
    "    retries = Retry(total=10, backoff_factor=0.1, status_forcelist=[500, 502, 503, 504])\n",
    "    if 'https' in url:\n",
    "        internet_status.mount('https://', HTTPAdapter(max_retries=retries))\n",
    "    elif 'http' in url:\n",
    "        internet_status.mount('http://', HTTPAdapter(max_retries=retries))\n",
    "    output_feedback = internet_status.get(url.format(Z=zoom_level, X=x, Y=y), verify=True)\n",
    "    if output_feedback.status_code == 200:\n",
    "        if loop_iterator == 0:\n",
    "            min_tile_index_width = x\n",
    "            min_tile_index_height = y\n",
    "\n",
    "        incoming_image = Image.open(BytesIO(output_feedback.content))\n",
    "        x = (x - min_tile_index_width) * 256\n",
    "        y = (y - min_tile_index_height) * 256\n",
    "        stitched_image.paste(im=incoming_image, box=(x, y))\n",
    "        incoming_image.close()\n",
    "        loop_iterator += 1\n",
    "        print(\"Number of downloaded images: \", loop_iterator, \" Percent Done: \", int((loop_iterator / total_number_of_tiles) * 100), \" %\", end=\"\\r\")\n",
    "\n",
    "# Here I will output the image to the desktop and not the browser. This is because the image can be several\n",
    "# - 100 megabytes big, depending on the given zoom level. Next, I save the image to the current working directory.\n",
    "# - Next, I load the image from memory into a numpy array, and close the image. I do this here becasue if the \n",
    "# - image is big then we can free up all of our memory for the numpy array. Next, I take the transform for our\n",
    "# - geo-referenced image. The transform operation will take the calcuated coordinates of out top left and bottom\n",
    "# - right tiles, and the width and height of the image. Next, I have to move the axis around so that our number\n",
    "# - of bands is in front of the tuple, rather than at the back. Since our image is RGB, the number of bands is 3,\n",
    "# - for three colors. Next, I take the image that is in the numpy array and open it with the rasterio call. Lastly,\n",
    "# - I write the image to disk, and print out where the file was saved to and the size of the image.\n",
    "stitched_image.show()\n",
    "file_path = os.path.join(os.getcwd(), \"Thy_Picture.tiff\")\n",
    "image_array = np.array(stitched_image)\n",
    "stitched_image.close()\n",
    "transform = rasterio.transform.from_bounds(first_long, second_lat, second_long, \n",
    "                                           first_lat, number_of_x_tiles * 256, number_of_y_tiles * 256)\n",
    "image_array = np.moveaxis(image_array, -1, 0)\n",
    "with rasterio.open(file_path, 'w', driver='GTiff', height=image_array.shape[1], \n",
    "                width=image_array.shape[2], count=image_array.shape[0], dtype=image_array.dtype,\n",
    "                crs='EPSG:3957', transform=transform) as dst:\n",
    "    dst.write(image_array, indexes=[1, 2, 3])\n",
    "\n",
    "print(file_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To add: \n",
    "    # Progress bar / pecentage\n",
    "    # Show the actual bounding box of what grabbed and then thee Tiles of what were grabbed."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}