cw_tiler_Demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Demo of CW-Tiler\n",
    "This demo will walk through the use of the tiler to create utm tiles with the SpaceNet Data Repository.\n",
    "\n",
    "We will be taking advantage of cloud optimized geotiffs.  For more information about SpaceNet visit https://spacenetchallenge.github.io/\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import base tools\n",
    "\n",
    "## Note, for mac osx compatability import something from shapely.geometry before importing fiona or geopandas\n",
    "## https://github.com/Toblerity/Shapely/issues/553  * Import shapely before rasterio or fioana\n",
    "from shapely import geometry\n",
    "import rasterio\n",
    "import random\n",
    "from cw_tiler import main\n",
    "from cw_tiler import utils\n",
    "from cw_tiler import vector_utils\n",
    "import numpy as np\n",
    "import os\n",
    "from tqdm import tqdm\n",
    "# Setting Certificate Location for Ubuntu/Mac OS locations (Rasterio looks for certs in centos locations)\n",
    "os.environ['CURL_CA_BUNDLE']='/etc/ssl/certs/ca-certificates.crt'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "## give location to SpaceNet 8-Band PanSharpened Geotiff on s3\n",
    "\n",
    "#spacenetPath = \"s3://spacenet-dataset/AOI_2_Vegas/srcData/rasterData/AOI_2_Vegas_MUL-PanSharpen_Cloud.tif\"\n",
    "spacenetPath = \"/home/dlindenbaum/datastorage/spacenet_sample/AOI_2_Vegas_MUL-PanSharpen_Cloud.tif\"\n",
    "osm_labels_path = \"./tests/fixtures/my-bucket/spacenet_test/las-vegas_nevada_osm_buildings.geojson\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Prep files for UTM\n",
    "with rasterio.open(spacenetPath) as src:\n",
    "\n",
    "    # Get Lat, Lon bounds of the Raster (src)\n",
    "    wgs_bounds = utils.get_wgs84_bounds(src)\n",
    "    \n",
    "    # Use Lat, Lon location of Image to get UTM Zone/ UTM projection\n",
    "    utm_crs = utils.calculate_UTM_crs(wgs_bounds)\n",
    "    \n",
    "    # Calculate Raster bounds in UTM coordinates \n",
    "    utm_bounds = utils.get_utm_bounds(src, utm_crs)\n",
    "\n",
    "## read vector file\n",
    "gdf = vector_utils.read_vector_file(osm_labels_path)\n",
    "gdf.head()\n",
    "gdf_utm = vector_utils.transformToUTM(gdf, utm_crs=utm_crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(652036.1651090909, 3998726.3249221705, 666254.6274175395, 4014808.742985228)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#print(gdf_utm.total_bounds)\n",
    "#utm_bounds\n",
    "\n",
    "geoBounds = geometry.box(*gdf_utm.total_bounds)\n",
    "rasterBounds = geometry.box(*utm_bounds)\n",
    "\n",
    "interBounds = geoBounds.intersection(rasterBounds)\n",
    "interBounds.bounds\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/dlindenbaum/miniconda3/envs/cw-tiler/lib/python3.6/site-packages/ipykernel/__main__.py:17: RasterioDeprecationWarning: Creation keywords stored on datasets by Rasterio versions < 1.0b1 will always be ignored in version 1.0. You may opt in to ignoring them now by setting RIO_IGNORE_CREATION_KWDS=TRUE in your environment.\n",
      "100%|██████████| 1476/1476 [1:33:10<00:00,  3.79s/it]\n"
     ]
    }
   ],
   "source": [
    "# open s3 Location\n",
    "\n",
    "# Each grid starting point will be spaced 400m apart\n",
    "stride_size_meters = 400\n",
    "\n",
    "# Each grid cell will be 400m on a side\n",
    "cell_size_meters   = 400\n",
    "\n",
    "# Specify the number of pixels in a tile cell_size_meters/tile_size_pixels == Pixel_Size_Meters\n",
    "tile_size_pixels   = 1200\n",
    "\n",
    "dataLocation = \"data\"\n",
    "imagePrefix = \"AOI_2_Vegas_MUL_{}_{}\"\n",
    "with rasterio.open(spacenetPath) as src:\n",
    "\n",
    "    \n",
    "    src_profile = src.profile\n",
    "    # Generate list of cells to read from utm_bounds \n",
    "    cells_list = main.calculate_analysis_grid(interBounds.bounds, stride_size_meters=stride_size_meters, cell_size_meters=cell_size_meters)\n",
    "\n",
    "    # select random cell\n",
    "    for cell_selection in tqdm(cells_list):\n",
    "        #random_cell = random.choice(cells_list)\n",
    "        ll_x, ll_y, ur_x, ur_y = cell_selection\n",
    "\n",
    "\n",
    "        # Get Tile from bounding box\n",
    "        tile, mask, window_transform = main.tile_utm(src, ll_x, ll_y, ur_x, ur_y, indexes=None, tilesize=tile_size_pixels, nodata=None, alpha=None,\n",
    "                     dst_crs=utm_crs)\n",
    "\n",
    "        #print(np.shape(tile))\n",
    "\n",
    "        ## Get Vector Information from bounding box\n",
    "        small_gdf = vector_utils.vector_tile_utm(gdf_utm, tile_bounds=[ll_x, ll_y, ur_x, ur_y])\n",
    "        if not small_gdf.empty:\n",
    "            small_gdf.to_file(os.path.join(dataLocation, imagePrefix.format(int(ll_x), int(ll_y))+\"_label.geojson\"), driver='GeoJSON')\n",
    "        \n",
    "        else:\n",
    "            open(os.path.join(dataLocation, imagePrefix.format(int(ll_x), int(ll_y))+\"_label.geojson\"), 'a').close()\n",
    "        \n",
    "        \n",
    "        img = vector_utils.rasterize_gdf(small_gdf,\n",
    "                                         src_shape=(tile_size_pixels, tile_size_pixels),\n",
    "                                         src_transform=window_transform,\n",
    "                                    )\n",
    "        \n",
    "        # update src_profile for writing new img\n",
    "        dst_profile = src_profile\n",
    "        dst_profile.update({'transform': window_transform,\n",
    "                    'crs': utm_crs,\n",
    "                    'width': tile_size_pixels,\n",
    "                    'height': tile_size_pixels,\n",
    "                            'count': 1,\n",
    "                            'dtype': rasterio.uint8\n",
    "\n",
    "                   })\n",
    "        ## write label to tiff\n",
    "        with rasterio.open(os.path.join(dataLocation, imagePrefix.format(int(ll_x), int(ll_y))+\"_label.tif\"),\n",
    "                           'w',\n",
    "                           **dst_profile) as dst:\n",
    "\n",
    "            dst.write(img, indexes=1)\n",
    "        \n",
    "        ## write image to geotiff\n",
    "        dst_profile = src_profile\n",
    "        dst_profile.update({'transform': window_transform,\n",
    "                    'crs': utm_crs,\n",
    "                    'width': tile_size_pixels,\n",
    "                    'height': tile_size_pixels,\n",
    "                            'count': 8,\n",
    "                            'dtype': rasterio.uint16\n",
    "                   })\n",
    "        \n",
    "        with rasterio.open(os.path.join(dataLocation, imagePrefix.format(int(ll_x), int(ll_y))+\"_image.tif\"), 'w',\n",
    "                           **dst_profile) as dst:\n",
    "\n",
    "            dst.write(tile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with concurrent.futures.ThreadPoolExecutor(\n",
    "                    max_workers=num_workers\n",
    "                ) as executor:\n",
    "\n",
    "                    # We map the compute() function over the raster\n",
    "                    # data generator, zip the resulting iterator with\n",
    "                    # the windows list, and as pairs come back we\n",
    "                    # write data to the destination dataset.\n",
    "                    for window, result in zip(\n",
    "                        windows, executor.map(compute, data_gen)\n",
    "                    ):\n",
    "                        dst.write(result, window=window)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "# Plot 1st Band of the Tile \n",
    "from  matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "#plt.imshow(img)\n",
    "plt.imshow(tile[1,:,:]) \n",
    "if not small_gdf.empty:\n",
    "    small_gdf.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "    ## Set for specific tile in Las Vegas\n",
    "    utm_crs = utils.calculate_UTM_crs([-115.30170, 36.15604, -115.30170, 36.15604])\n",
    "    gdf = vector_utils.read_vector_file(geojsonPath)\n",
    "    gdf.head()\n",
    "    gdf = vector_utils.transformToUTM(gdf, utm_crs=utm_crs)\n",
    "\n",
    "    utmX, utmY = 658029, 4006947\n",
    "    ll_x = utmX\n",
    "    ur_x = utmX + 500\n",
    "    ll_y = utmY\n",
    "    ur_y = utmY + 500\n",
    "    stride_size_meters = 300\n",
    "    cell_size_meters = 400\n",
    "    tile_size_pixels = 1600\n",
    "\n",
    "    address = spacenetPath\n",
    "\n",
    "    with rasterio.open(address) as src:\n",
    "\n",
    "        src_profile = src.profile\n",
    "\n",
    "\n",
    "        tile, mask, window_transform = main.tile_utm(src, ll_x, ll_y, ur_x, ur_y,\n",
    "                                                     indexes=None,\n",
    "                                                     tilesize=tile_size_pixels,\n",
    "                                                     nodata=None,\n",
    "                                                     alpha=None,\n",
    "                                                     dst_crs=utm_crs\n",
    "                                                     )\n",
    "\n",
    "        print(np.shape(tile))\n",
    "        assert np.shape(tile) == (8, tile_size_pixels, tile_size_pixels)\n",
    "\n",
    "        small_gdf = vector_utils.vector_tile_utm(gdf, tile_bounds=[ll_x, ll_y, ur_x, ur_y])\n",
    "        print(small_gdf.shape)\n",
    "        small_gdf.to_file(os.path.join(PREFIX, \"testTiff_Label.geojson\"), driver='GeoJSON')\n",
    "        img = vector_utils.rasterize_gdf(small_gdf,\n",
    "                                         src_shape=(tile_size_pixels, tile_size_pixels),\n",
    "                                         src_transform=window_transform,\n",
    "                                    )\n",
    "        print(\"Label Count Burn {}:\".format(np.sum(img)))\n",
    "        assert img.shape == (tile_size_pixels, tile_size_pixels)\n",
    "\n",
    "\n",
    "        dst_profile = src_profile\n",
    "        dst_profile.update({'transform': window_transform,\n",
    "                    'crs': utm_crs,\n",
    "                    'width': tile_size_pixels,\n",
    "                    'height': tile_size_pixels,\n",
    "                            'count': 1,\n",
    "                            'dtype': rasterio.uint8\n",
    "\n",
    "                   })\n",
    "\n",
    "        with rasterio.open(os.path.join(PREFIX, \"testTiff_Label.tif\"), 'w',\n",
    "                                        **dst_profile) as dst:\n",
    "\n",
    "            dst.write(img, indexes=1)\n",
    "\n",
    "        dst_profile = src_profile\n",
    "        dst_profile.update({'transform': window_transform,\n",
    "                    'crs': utm_crs,\n",
    "                    'width': tile_size_pixels,\n",
    "                    'height': tile_size_pixels,\n",
    "                            'count': 8,\n",
    "                            'dtype': rasterio.uint16\n",
    "                   })\n",
    "\n",
    "        with rasterio.open(os.path.join(PREFIX, \"testTiff_Image.tif\"), 'w',\n",
    "                           **dst_profile) as dst:\n",
    "\n",
    "            dst.write(tile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas as gpd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "small_gdf = gpd.GeoDataFrame(geometry=[])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "small_gdf.area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:cw-tiler]",
   "language": "python",
   "name": "conda-env-cw-tiler-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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}