Height Extractor Code.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "import os\n",
    "import pandas as pd\n",
    "import geopandas as gpd #important\n",
    "import rasterio\n",
    "import matplotlib.pyplot as plt\n",
    "import operator\n",
    "import seaborn as sns\n",
    "import psycopg2\n",
    "import shapely \n",
    "import laspy #las open\n",
    "from laspy.file import File as read_las #open las \n",
    "from shapely.geometry import Point #convert las to gpd\n",
    "import descartes #map figures\n",
    "import rtree #spatial merge\n",
    "import statistics #median"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#1. READ LAS#\n",
    "in_las = read_las('lassample.las', mode='r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#2. CONVERTING LAS TO PANDAS#\n",
    "\n",
    "#Import LAS into numpy array (X=raw integer value; x=scaled float value)\n",
    "lidar_points = np.array((in_las.x,in_las.y,in_las.z,in_las.intensity)).transpose()\n",
    "\n",
    "#Transform to pandas DataFrame\n",
    "lidar_df=pd.DataFrame(lidar_points, columns = ['X coordinate', 'Y coordinate', 'Z coordinate', 'Intensity'])\n",
    "\n",
    "#Transform to geopandas GeoDataFrame\n",
    "Geometry = [Point(xy) for xy in zip(in_las.x,in_las.y)]\n",
    "lidar = gpd.GeoDataFrame(lidar_df, crs = {'init' :'epsg:32643'}, geometry=Geometry) #set correct spatial reference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#3. READ SHAPEFILE#\n",
    "shapefile = gpd.read_file(\"shapefilesample.shp\")\n",
    "\n",
    "#make id the index\n",
    "shapefile = shapefile.set_index(shapefile.id)\n",
    "shapefile  = shapefile.drop(['id'], axis = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#4. SPATIAL MERGE LIDAR AND SHAPEFILE BY GEOMETRY COLUMNS# \n",
    "#Keep building id (stored as index_right)\n",
    "buildings = gpd.sjoin(lidar, shapefile) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#5. CALCULATE MEDIAN HEIGHT OF BUILDING BY GROUPING BY BUILDING ID# \n",
    "\n",
    "building_height_input = buildings.groupby('index_right')['Z coordinate'].median() #add filters here if needed\n",
    "building_height = pd.DataFrame(building_height_input, index = None )\n",
    "building_height = building_height.rename_axis(index='ID').rename(columns={'Z coordinate':'Roof Height'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#6. CREATE BUFFER FROM THE BUILDING EDGES#\n",
    "\n",
    "#create 1 metre buffer\n",
    "buffer_line1a = shapefile['geometry'].buffer(distance = 2) #returns Geoseries\n",
    "buffer_line1 = gpd.GeoDataFrame(buffer_line1a, columns = ['geometry'], crs = {'init' :'epsg:32643'}) #convert to Geodatabase\n",
    "\n",
    "#create 2 metre buffer\n",
    "buffer_line2a = shapefile['geometry'].buffer(distance = 1)\n",
    "buffer_line2 = gpd.GeoDataFrame(buffer_line2a, columns = ['geometry'], crs = {'init' :'epsg:32643'})\n",
    "\n",
    "#Take difference between two buffers\n",
    "buffer_line = buffer_line1['geometry'].difference(buffer_line2)\n",
    "buffer_line = gpd.GeoDataFrame(buffer_line, columns = ['geometry'], crs = {'init' :'epsg:32643'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#7. SPATIAL MERGE LIDAR AND BUFFER_LINE BY GEOMETRY COLUMNS# \n",
    "#Keep building id (stored as index_right)\n",
    "border = gpd.sjoin(lidar, buffer_line)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#8. CALCULATE GROUND LEVEL OF LIDAR Z POINTS IN BORDER BY GROUPING BY BUILDING ID #\n",
    "#(version with 0.025 quantile taken as the parameter)\n",
    "\n",
    "ground_height_input = border.groupby('index_right')['Z coordinate'].quantile(q=0.025)\n",
    "#if (ground_height_input)\n",
    "ground_height = pd.DataFrame(ground_height_input)\n",
    "ground_height = ground_height.rename_axis(index='ID').rename(columns={'Z coordinate':'Ground Level'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#8. CALCULATE GROUND LEVEL OF LIDAR Z POINTS IN BORDER BY GROUPING BY BUILDING ID#\n",
    "#(version with minimum taken as the parameter)\n",
    "\n",
    "ground_height_input2 = border.groupby('index_right')['Z coordinate'].min()\n",
    "ground_height2 = pd.DataFrame(ground_height_input2)\n",
    "ground_height2 = ground_height2.rename_axis(index='ID').rename(columns={'Z coordinate':'Ground Level'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#9. CREATE DATAFRAME OF BUILDING ID AND BUILDING HEIGHT#\n",
    "\n",
    "height_difference1 = building_height['Roof Height'] - ground_height['Ground Level']\n",
    "#height_difference.dropna(how = 'any', inplace = True) #drops null values \n",
    "height_difference = pd.DataFrame({'Building Height' : height_difference1})\n",
    "height_difference = height_difference.rename_axis(index='ID')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "#10. MERGE INTO FINAL DETAILED TABLE#\n",
    "\n",
    "out_final = shapefile.join(height_difference)\n",
    "out_final = out_final.dropna(how = 'any', inplace = False)\n",
    "out_final['Building ID'] = out_final.index\n",
    "out_final['Roof Height'] = building_height['Roof Height']\n",
    "out_final['Ground Level'] = ground_height['Ground Level']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "#11. CSV EXTRACTION#\n",
    "\n",
    "out_final.to_csv('height.csv',index = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {},
   "outputs": [],
   "source": [
    "#11. SHP EXTRACTION#\n",
    "\n",
    "out_final.to_file(\"height_out_with_details.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "#analysis of lidar points collected on building tops#\n",
    "building_summary = buildings.groupby('index_right')['Z coordinate'].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "#analysis of lidar points collected on ground buffer zone#\n",
    "ground_summary = border.groupby('index_right')['Z coordinate'].describe()"
   ]
  }
 ],
 "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}