Example_Code_Template 3D_Houses_Belgium.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot 3D houses in Belgium"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Address"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# enter an address in Belgium with the following format \n",
    "# \"road_name number postcode name_of_the_city\"\n",
    "# example address = \"Molenstraat 129 1570 Galmaarden\"\n",
    "\n",
    "address = input(\"Enter an address in Belgium:\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pandas as pd\n",
    "\n",
    "from glob import glob  \n",
    "from natsort import natsorted, ns\n",
    "\n",
    "import geopy\n",
    "import folium\n",
    "from geopy.geocoders import Nominatim\n",
    "from geopy.extra.rate_limiter import RateLimiter\n",
    "\n",
    "import rasterio  \n",
    "from rasterio.plot import show\n",
    "from pyproj import Transformer\n",
    "from rasterio.windows import Window\n",
    "import rioxarray \n",
    "\n",
    "import imageio\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "from osgeo import gdal\n",
    "from mayavi import mlab\n",
    "\n",
    "import webbrowser\n",
    "from IPython.display import Image\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# setting your dataset path - Flandre dataset\n",
    "# flandre_path = '/Users/adamtky/Desktop/BeCode/3. Project/Project 3 - 3D Houses'\n",
    "\n",
    "flandre_path = '/your_dataset_path/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# setting your dataset path - Wallonie dataset\n",
    "# wallonie_path = '/Volumes/Seagate Expansion Drive/Project 3 - 3D Houses_Wallonie'\n",
    "\n",
    "wallonie_path = '/your_dataset_path/'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Search all tifs and sort files on laptop ( Flandre & Brussels )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def tif_file(path):\n",
    "    \n",
    "    tif_file =[]\n",
    "    \n",
    "    # using glob library to get all the file with .tif\n",
    "    files = glob(path,recursive = True) \n",
    "    for file in files: \n",
    "        tif_file.append(file)\n",
    "        \n",
    "    # sort files with number in the file (ascending order) using natsorted library\n",
    "    tif_file = natsorted(tif_file, alg=ns.IGNORECASE)     \n",
    "    return tif_file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path = './DSM/**/*.tif'\n",
    "DSM_tif = tif_file(path)\n",
    "DSM_tif[:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path = './DTM/**/*.tif'\n",
    "DTM_tif = tif_file(path)\n",
    "DTM_tif[:3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Search all tifs and sort files on hard drive ( Wallanie )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(wallonie_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path_wal = './DSM/*.tif'\n",
    "DSM_tif_wal = tif_file(path_wal)\n",
    "DSM_tif_wal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path_wal = './DTM/*.tif'\n",
    "DTM_tif_wal = tif_file(path_wal)\n",
    "DTM_tif_wal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "### Geocoding - single address"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# to get the longtitude and latitude of the address entered & plot address on a map\n",
    "\n",
    "def lat_long(address):\n",
    "    \n",
    "    # to get the longtitude and latitude \n",
    "    \n",
    "    locator = Nominatim(user_agent=\"myGeocoder\")\n",
    "    location = locator.geocode(address)\n",
    "    house_lat_long = [location.latitude, location.longitude]\n",
    "    return house_lat_long\n",
    "\n",
    "def house_map(func):\n",
    "    \n",
    "    # to plot address \n",
    "    \n",
    "    house_map = folium.Map(location=func,zoom_start=150)\n",
    "    folium.Marker(location=func,popup=lat_long(address)).add_to(house_map)\n",
    "    return house_map \n",
    "\n",
    "house_map(lat_long(address))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convert longtitude & latitude "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# longtitude and latitude of the address entered\n",
    "lat,lon = lat_long(address)\n",
    "lat,lon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def transform(lon, lat):\n",
    "    \n",
    "    # transform to Belgium 'EPSG:31370' coordinate\n",
    "    transformer = Transformer.from_crs(\"EPSG:4326\", crs_to = 'EPSG:31370' ,always_xy=True)\n",
    "    xx, yy = transformer.transform(lon, lat)\n",
    "    print(f'longtitude:{xx} latitude:{yy}' )\n",
    "    \n",
    "    return xx,yy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# coordinate EPSG:31370\n",
    "xx,yy = transform(lon, lat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using CSV to locate coordinate EPSG:31370"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# enter address_path\n",
    "# i.e address_path = '/Users/adamtky/Desktop/BeCode/3. Project/Project 3 - 3D Houses/Data/Belgium_Address.csv'\n",
    "\n",
    "address_path ='your_path/Belgium_Address.csv'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_coordinate(path,address):\n",
    "    \n",
    "    index = []\n",
    "    \n",
    "    df = pd.read_csv(path)\n",
    "    \n",
    "    for i in range(df.shape[0]):\n",
    "\n",
    "        if df['Address'][i] == address :\n",
    "            index.append(i)\n",
    "            break\n",
    "        continue\n",
    "\n",
    "    xx, yy = df['EPSG:31370_x'][index][index[0]] ,df['EPSG:31370_y'][index][index[0]] \n",
    "\n",
    "    lat, lon = df['EPSG:4326_lat'][index][index[0]] ,df['EPSG:4326_lon'][index][index[0]] \n",
    "    \n",
    "    return (xx, yy) , (lat, lon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(xx, yy) , (lat, lon) = get_coordinate(address_path,address)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Locate tif for given longtitude & latitude"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(flandre_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create all bounding box from tifs\n",
    "\n",
    "def bounds(tifs):\n",
    "    \n",
    "    bounds = []\n",
    "    \n",
    "    for i in tifs:\n",
    "        src = rasterio.open(i)\n",
    "        bounds.append(src.bounds)\n",
    "        \n",
    "    return  bounds "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "src_bounds_flb = bounds(DSM_tif)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "src_bounds_flb[:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(wallonie_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "src_bounds_wal = bounds(DSM_tif_wal)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "src_bounds_wal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check which tif contains the location in DSM\n",
    "\n",
    "def check_DSM_flb(xx,yy):\n",
    "    \n",
    "    scr_path = []\n",
    "\n",
    "    for i,bound in enumerate(src_bounds_flb,1):\n",
    "\n",
    "        if (xx >= bound[0] and xx <= bound[2]) & \\\n",
    "           (yy >= bound[1] and yy <= bound[3]):\n",
    "\n",
    "            scr_path.append('./DSM/DHMVIIDSMRAS1m_k'+ str(i) +'/GeoTIFF/DHMVIIDSMRAS1m_k'+ str(i) + '.tif')  \n",
    "\n",
    "            print('The house is in this tif :', 'DHMVIIDSMRAS1m_k' + str(i) + '.tif')\n",
    "        else:\n",
    "            None\n",
    "            \n",
    "    return scr_path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dsm_path_flb = check_DSM_flb(xx,yy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dsm_path_flb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check which tif contains the location in DSM\n",
    "\n",
    "def check_DSM_wal(xx,yy):\n",
    "\n",
    "    scr_path = []\n",
    "\n",
    "    for i,bound in enumerate(src_bounds_wal,1):\n",
    "\n",
    "        if (xx >= bound[0] and xx <= bound[2]) & \\\n",
    "           (yy >= bound[1] and yy <= bound[3]):\n",
    "\n",
    "            scr_path.append('/Volumes/Seagate Expansion Drive/Project 3 - 3D Houses_Wallonie/DSM/dsm_'+ str(i) + '.tif')  \n",
    "\n",
    "            print('The house is in this tif :', '/Volumes/Seagate Expansion Drive/Project 3 - 3D Houses_Wallonie/DSM/dsm_' + str(i) + '.tif')\n",
    "            \n",
    "    return scr_path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dsm_path_wal = check_DSM_wal(xx,yy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dsm_path_wal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check which tif contains the location in DTM\n",
    "\n",
    "\n",
    "def check_DTM_flb(xx,yy):\n",
    "    \n",
    "    scr_path = []\n",
    "\n",
    "    for i,bound in enumerate(src_bounds_flb,1):\n",
    "\n",
    "        if (xx >= bound[0] and xx <= bound[2]) & \\\n",
    "           (yy >= bound[1] and yy <= bound[3]):\n",
    "\n",
    "            scr_path.append('./DTM/DHMVIIDTMRAS1m_k'+ str(i) +'/GeoTIFF/DHMVIIDTMRAS1m_k'+ str(i) + '.tif')  \n",
    "\n",
    "            print('The house is in this tif :', 'DHMVIIDTMRAS1m_k' + str(i) + '.tif')\n",
    "        else:\n",
    "            None\n",
    "            \n",
    "    return scr_path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dtm_path_flb = check_DTM_flb(xx,yy)\n",
    "dtm_path_flb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check which tif contains the location in DTM\n",
    "\n",
    "def check_DTM_wal(xx,yy):\n",
    "\n",
    "    scr_path = []\n",
    "\n",
    "    for i,bound in enumerate(src_bounds_wal,1):\n",
    "\n",
    "        if (xx >= bound[0] and xx <= bound[2]) & \\\n",
    "           (yy >= bound[1] and yy <= bound[3]):\n",
    "\n",
    "            scr_path.append('/Volumes/Seagate Expansion Drive/Project 3 - 3D Houses_Wallonie/DTM/dtm_'+ str(i) + '.tif')  \n",
    "\n",
    "            print('The house is in this tif :', '/Volumes/Seagate Expansion Drive/Project 3 - 3D Houses_Wallonie/DTM/dtm_' + str(i) + '.tif')\n",
    "            \n",
    "    return scr_path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dtm_path_wal = check_DTM_wal(xx,yy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clip the location of the house from tif"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dsm_path = []\n",
    "\n",
    "if not dsm_path_flb:\n",
    "    dsm_path.append(dsm_path_wal[0])\n",
    "    os.chdir(wallonie_path)\n",
    "else:\n",
    "    dsm_path.append(dsm_path_flb[0])\n",
    "    os.chdir(flandre_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dtm_path = []\n",
    "\n",
    "if not dtm_path_flb:\n",
    "    dtm_path.append(dtm_path_wal[0])\n",
    "    os.chdir(wallonie_path)\n",
    "else:\n",
    "    dtm_path.append(dtm_path_flb[0])\n",
    "    os.chdir(flandre_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dsm_path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dtm_path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def clip_dsm(path,window_size:int):\n",
    " \n",
    "    xds = rioxarray.open_rasterio(path,masked=True,chunks=True)\n",
    "    \n",
    "    # set window size\n",
    "    n = window_size\n",
    "\n",
    "    # create coordinates and geometries\n",
    "    coor1,coor2 = [(xx-n),(yy+n)],[(xx+n),(yy+n)]\n",
    "\n",
    "    coor3,coor4 = [(xx+n),(yy-n)] ,[(xx-n),(yy-n)]\n",
    "\n",
    "    geometries = [ {'type': 'Polygon', 'coordinates': [[coor1,coor2, coor3,coor4,coor1 ]]}]\n",
    "\n",
    "    # clip the image as per the geometries size\n",
    "    clipped = xds.rio.clip(geometries)\n",
    "    \n",
    "    # change directory\n",
    "    os.chdir(flandre_path)\n",
    "    \n",
    "    # save clip\n",
    "    clip_dsm = clipped.rio.to_raster(address+\"_clipped_dsm.tif\",tiled=True, dtype=\"int32\")\n",
    "\n",
    "    return clipped.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def clip_dtm(path,window_size:int):\n",
    " \n",
    "    xds = rioxarray.open_rasterio(path,masked=True,chunks=True)\n",
    "    \n",
    "    # set window size\n",
    "    n = window_size\n",
    "\n",
    "    # create coordinates and geometries\n",
    "    coor1,coor2 = [(xx-n),(yy+n)],[(xx+n),(yy+n)]\n",
    "\n",
    "    coor3,coor4 = [(xx+n),(yy-n)] ,[(xx-n),(yy-n)]\n",
    "\n",
    "    geometries = [ {'type': 'Polygon', 'coordinates': [[coor1,coor2, coor3,coor4,coor1 ]]}]\n",
    "\n",
    "    # clip the image as per the geometries size\n",
    "    clipped = xds.rio.clip(geometries)\n",
    "\n",
    "    # change directory\n",
    "    os.chdir(flandre_path)\n",
    "    \n",
    "    # save clip\n",
    "    clip_dtm = clipped.rio.to_raster(address +\"_clipped_dtm.tif\",tiled=True, dtype=\"int32\")\n",
    "\n",
    "    return clipped.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "clip_dsm(dsm_path[0],20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clip_dtm(dtm_path[0],20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def chm_tif():\n",
    "    \n",
    "    # open the digital terrain model \n",
    "    with rasterio.open(address +'_clipped_dtm.tif') as src:\n",
    "        lidar_dtm_im = src.read(1, masked=True)\n",
    "        dtm_meta = src.profile\n",
    "        \n",
    "    # open the digital surface model \n",
    "    with rasterio.open(address +'_clipped_dsm.tif') as src:\n",
    "        lidar_dsm_im = src.read(1, masked=True)\n",
    "        dsm_meta = src.profile\n",
    "    \n",
    "    # calculate canopy height model\n",
    "    lidar_chm = lidar_dsm_im - lidar_dtm_im\n",
    "    \n",
    "    # change directory\n",
    "    os.chdir(flandre_path)\n",
    "    \n",
    "    # save chm clipped \n",
    "    with rasterio.open(address +'_clipped_chm.tif', 'w', **dsm_meta) as ff:\n",
    "        ff.write(lidar_chm,1)\n",
    "        \n",
    "    chm_tif = address +'_clipped_chm.tif'\n",
    "    \n",
    "    return chm_tif"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create 3D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def House_3D(tif,size:int,azim=215):\n",
    "\n",
    "    chm = imageio.imread(tif)\n",
    "\n",
    "    nx,ny = chm.shape\n",
    "\n",
    "    x = np.linspace(0, size*2, nx)\n",
    "    y = np.linspace(0, size*2, ny)\n",
    "\n",
    "    yv,xv = np.meshgrid(x, y)\n",
    "\n",
    "    fig = plt.figure(figsize=(15,15))\n",
    "    ax = fig.add_subplot(111, projection='3d')\n",
    "\n",
    "    chm3d=ax.plot_surface(xv,yv,chm,cmap='plasma',linewidth=0)\n",
    "    ax.set_title('CHM(Canopy Height Model)')\n",
    "    ax.set_xlabel('Distance (m)')\n",
    "    ax.set_ylabel('Distance (m)')\n",
    "    ax.set_zlabel('Elevation (m)')\n",
    "    ax.view_init(azim=azim)\n",
    "        \n",
    "    fig.colorbar(chm3d, shrink=0.3, aspect=10);\n",
    "    fig.savefig(address +'_3D.png', dpi=200) \n",
    "    \n",
    "    return  plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "House_3D(chm_tif(),20,45)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# to view actual house on google map\n",
    "url = 'https://www.google.com.my/maps/place/'+str(lat)+','+str(lon)\n",
    "webbrowser.open(url)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting the CHM with Mayavi2\n",
    "\n",
    "ds = gdal.Open(chm_tif())\n",
    "data = ds.ReadAsArray()\n",
    "\n",
    "data = data.astype(np.float32)\n",
    "\n",
    "mlab.figure(size=(640, 800), bgcolor=(0.16, 0.28, 0.46))\n",
    "\n",
    "surf = mlab.surf(data)\n",
    "\n",
    "mlab.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compare google map image of the house\n",
    "\n",
    "# Image(address+\".png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mayavi image\n",
    "# Image(address+\"_mayavi.png\")"
   ]
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}