lasAnalysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "630c3eeb-36aa-4a77-89cd-575c8a58ab94",
   "metadata": {},
   "source": [
    "To run this notebook you need to have several packages installed. One way to install in through Anaconda with the following:\n",
    "\n",
    "`conda env create -n jack -f environment.yml`\n",
    "\n",
    "`conda activate jack`\n",
    "\n",
    "You can work with the las data two ways:\n",
    "1. Use the las2rec function to read the data into a recArray.\n",
    "2. Access the point directly using laspy, such as las.x or las.y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2d3cad9d-7400-4947-b073-29e80876ab0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import laspy\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from pyproj import crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "9ece75ef-3680-4f34-b13e-cd7152a23f3e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6645595\n",
      "4904697\n"
     ]
    }
   ],
   "source": [
    "def las2rec(infile):\n",
    "    \"\"\"\n",
    "    Function to read in a laz file and return as a recArray.\n",
    "    \"\"\"\n",
    "    with laspy.open(infile) as fh:\n",
    "        las = fh.read()\n",
    "    las_data = np.rec.fromarrays([las.return_num, las.num_returns, las.gps_time, las.intensity,las.classification, las.x, las.y, las.z],\n",
    "                                 names=['RETURN_NUMBER', 'NUMBER_OF_RETURNS', 'TIMESTAMP', 'INTENSITY', 'CLASSIFICATION', 'X', 'Y', 'Z'],\n",
    "                                 formats=['u1', 'u1', '<f8', '<i4', 'u1', '<f8', '<f8', '<f8'])\n",
    "    \n",
    "    las_data.sort(order='TIMESTAMP')\n",
    "    \n",
    "    print(las_data['X'].size)\n",
    "    \n",
    "    return las_data\n",
    "\n",
    "# Read data from las file into recArray \n",
    "lasFile1 = r'E:\\lid1_lid2_comparison\\input_data\\baradine_201407_lid1\\laz\\apl1dr_x692000ys6570000_20140703_ba1m5_pBARAL1.laz'\n",
    "lasData1 = las2rec(lasFile1)\n",
    "\n",
    "lasFile2 = r'E:\\lid1_lid2_comparison\\reduced_density_data\\baradine_201407_lid1\\laz_131\\apl1dr_x692000ys6570000_20140703_ba1m5_pBARAL1.laz'\n",
    "lasData2 = las2rec(lasFile2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "864885e5-1070-4a0c-97ef-c69c67b9faa0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Original data\n",
      "min and max X: 692000.000000 693999.999000\n",
      "min and max Y: 6568000.000000 6569999.999000\n",
      "min and max Z: 294.694000 349.857000\n",
      "Filtered data\n",
      "min and max X: 692000.000000 693999.999000\n",
      "min and max Y: 6568000.000000 6569999.999000\n",
      "min and max Z: 294.694000 349.857000\n"
     ]
    }
   ],
   "source": [
    "# Print out the bounding box\n",
    "print(\"Original data\")\n",
    "print('min and max X: %f %f'%(np.min(lasData1['X']), np.max(lasData1['X'])))\n",
    "print('min and max Y: %f %f'%(np.min(lasData1['Y']), np.max(lasData1['Y'])))\n",
    "print('min and max Z: %f %f'%(np.min(lasData1['Z']), np.max(lasData1['Z'])))\n",
    "print(\"Filtered data\")\n",
    "print('min and max X: %f %f'%(np.min(lasData1['X']), np.max(lasData1['X'])))\n",
    "print('min and max Y: %f %f'%(np.min(lasData1['Y']), np.max(lasData1['Y'])))\n",
    "print('min and max Z: %f %f'%(np.min(lasData1['Z']), np.max(lasData1['Z'])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "cc92b78b-94a7-45a4-b23b-4575a4f810e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotReturnNumberCount(las_data):\n",
    "    \"\"\"\n",
    "    Function to make a histogram of return numbers from a las recArray.\n",
    "    \"\"\"\n",
    "    returnNumbers = las_data['RETURN_NUMBER']\n",
    "    returnNumberList = np.unique(returnNumbers)\n",
    "    count = []\n",
    "    for n in returnNumberList:\n",
    "        count.append(np.sum(returnNumbers == n))\n",
    "        print(n, np.sum(returnNumbers == n))\n",
    "    fig = plt.figure(1)\n",
    "    fig.set_size_inches((8, 2))\n",
    "    ax = plt.axes([0.2, 0.75, 0.2, 0.75])  \n",
    "    ax.bar(returnNumberList, count)\n",
    "    ax.set_xlabel('Return numbers')\n",
    "    ax.set_ylabel('Count')\n",
    "    ax.set_xticks(returnNumberList)\n",
    "    ax.ticklabel_format(style='plain')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "04a32c27-aeb4-4bac-a41c-0650b80d47a6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 6329403\n",
      "2 287307\n",
      "3 28371\n",
      "4 514\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQ8AAADVCAYAAABJ/XXDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAi9UlEQVR4nO3df1RUdf7H8efIjwEJRgQBx1jB1iUJ3Qw8hubiT7Akc2vLPRQrm7GaCCK6GfZ1/dEuWGvqKkbZD91+GO2u0bJHJJAKc4VEgnRc1DxBYoLYNg5oBAj3+4eHu47jzys0Y74f53DOzue+59733LXX3Llz72d0iqIoCCHENepl7waEEDcmCQ8hhCYSHkIITSQ8hBCaSHgIITSR8BBCaCLhIYTQRMJDCKGJs70bENY6Ozs5fvw4np6e6HQ6e7cjbjKKotDc3IzRaKRXr8sfW0h4OJjjx48TGBho7zbETa6uro5bb731sjUSHg7G09MTOPd/npeXl527ETebpqYmAgMD1X+HlyPh4WC6Pqp4eXlJeAi7uZqPzHLCVAihiYSHEEITCQ8hhCYSHkIITeSE6Q0o6Oltdtlu7copdtmucExy5CGE0ETCQwihid3D4+uvv+axxx7Dx8eH3r17c+edd1JRUaEuVxSFZcuWYTQacXd3Z+zYsRw4cMBqHa2trSQnJ+Pr64uHhwdTp07l2LFjVjVms5n4+HgMBgMGg4H4+HhOnTplVXP06FHuv/9+PDw88PX1JSUlhba2Nqua/fv3ExUVhbu7OwMGDGDFihVcOA1sSUkJ4eHhuLm5MWjQIF566aVu2FNCOBa7hofZbGb06NG4uLiwfft2/vOf//DCCy/Qp08fteb5559n9erVZGVlUV5eTkBAAJMmTaK5uVmtSU1NJTc3l5ycHHbt2sXp06eJjY2lo6NDrYmLi6OqqoqCggIKCgqoqqoiPj5eXd7R0cGUKVM4c+YMu3btIicnh61bt7JgwQK1pqmpiUmTJmE0GikvL2f9+vWsWrWK1atXqzU1NTXcd999jBkzhsrKShYvXkxKSgpbt27tob0ohH3o7Dl7+tNPP82///1vPvnkk4suVxQFo9FIamoqixYtAs4dZfj7+/Pcc88xa9YsLBYL/fr1480332T69OnA/+4Pyc/PJyYmhurqakJDQykrK2PkyJEAlJWVERkZycGDBwkJCWH79u3ExsZSV1eH0WgEICcnh4SEBBobG/Hy8iI7O5v09HROnDiBXq8HYOXKlaxfv55jx46h0+lYtGgReXl5VFdXq69j9uzZfP7555SWll5xnzQ1NWEwGLBYLJe8wlROmIqecjX//rrY9cgjLy+PiIgIHn74Yfz8/Bg+fDivvPKKurympoaGhgaio6PVMb1eT1RUFLt37wagoqKC9vZ2qxqj0UhYWJhaU1paisFgUIMD4O6778ZgMFjVhIWFqcEBEBMTQ2trq/oxqrS0lKioKDU4umqOHz9ObW2tWnN+L101e/fupb293WYftLa20tTUZPUnxI3AruHx5Zdfkp2dzeDBg/nggw+YPXs2KSkpvPHGGwA0NDQA4O/vb/U8f39/dVlDQwOurq54e3tftsbPz89m+35+flY1F27H29sbV1fXy9Z0Pb5SzdmzZ/nmm29sesjMzFTPwxgMBrmjVtww7BoenZ2d3HXXXWRkZDB8+HBmzZpFYmIi2dnZVnUX3qSjKMoVb9y5sOZi9d1R0/Wp71pruqSnp2OxWNS/urq6y74uIRyFXcOjf//+hIaGWo0NGTKEo0ePAhAQEAD87129S2Njo/ruHhAQQFtbG2az+bI1J06csNn+yZMnrWou3I7ZbKa9vf2yNY2NjQBXrHF2dsbHx8emB71er95BK3fSihuJXcNj9OjRHDp0yGrs8OHDDBw4EIDg4GACAgIoKipSl7e1tVFSUsKoUaMACA8Px8XFxaqmvr4ek8mk1kRGRmKxWNizZ49a8+mnn2KxWKxqTCYT9fX1ak1hYSF6vZ7w8HC1ZufOnVZf3xYWFmI0GgkKClJrzu+lqyYiIgIXFxdtO0oIB2TX8Jg/fz5lZWVkZGRw5MgRtmzZwsaNG0lKSgLOHeanpqaSkZFBbm4uJpOJhIQEevfuTVxcHAAGg4GZM2eyYMECiouLqays5LHHHmPo0KFMnDgROHc0M3nyZBITEykrK6OsrIzExERiY2MJCQkBIDo6mtDQUOLj46msrKS4uJiFCxeSmJioHg3ExcWh1+tJSEjAZDKRm5tLRkYGaWlp6keS2bNn89VXX5GWlkZ1dTWvv/46r732GgsXLvyhd68QPcqu97aMGDGC3Nxc0tPTWbFiBcHBwaxdu5ZHH31UrXnqqadoaWlhzpw5mM1mRo4cSWFhodVMR2vWrMHZ2ZlHHnmElpYWJkyYwObNm3FyclJr3n77bVJSUtRvQqZOnUpWVpa63MnJiW3btjFnzhxGjx6Nu7s7cXFxrFq1Sq0xGAwUFRWRlJREREQE3t7epKWlkZaWptYEBweTn5/P/Pnz2bBhA0ajkXXr1vHQQw/1yD4Uwl7sep2HsCXXeQh7umGu8xBC3LgkPIQQmkh4CCE0kfAQQmgi4SGE0ETCQwihiYSHEEITCQ8hhCYSHkIITSQ8hBCaSHgIITSR8BBCaCLhIYTQRMJDCKGJhIcQQhMJDyGEJhIeQghNJDyEEJpIeAghNJHwEEJoIuEhhNBEwkMIoYmEhxBCEwkPIYQmEh5CCE0kPIQQmkh4CCE0kfAQQmgi4SGE0ETCQwihiYSHEEITCQ8hhCYSHkIITSQ8hBCaOEx4ZGZmotPpSE1NVccURWHZsmUYjUbc3d0ZO3YsBw4csHpea2srycnJ+Pr64uHhwdSpUzl27JhVjdlsJj4+HoPBgMFgID4+nlOnTlnVHD16lPvvvx8PDw98fX1JSUmhra3Nqmb//v1ERUXh7u7OgAEDWLFiBYqiWNWUlJQQHh6Om5sbgwYN4qWXXrr+nSOEA3KI8CgvL2fjxo0MGzbMavz5559n9erVZGVlUV5eTkBAAJMmTaK5uVmtSU1NJTc3l5ycHHbt2sXp06eJjY2lo6NDrYmLi6OqqoqCggIKCgqoqqoiPj5eXd7R0cGUKVM4c+YMu3btIicnh61bt7JgwQK1pqmpiUmTJmE0GikvL2f9+vWsWrWK1atXqzU1NTXcd999jBkzhsrKShYvXkxKSgpbt27tid0mhF3plAvfOq/CoEGDKC8vx8fHx2r81KlT3HXXXXz55ZdXva7Tp09z11138eKLL/LHP/6RO++8k7Vr16IoCkajkdTUVBYtWgScO8rw9/fnueeeY9asWVgsFvr168ebb77J9OnTATh+/DiBgYHk5+cTExNDdXU1oaGhlJWVMXLkSADKysqIjIzk4MGDhISEsH37dmJjY6mrq8NoNAKQk5NDQkICjY2NeHl5kZ2dTXp6OidOnECv1wOwcuVK1q9fz7Fjx9DpdCxatIi8vDyqq6vV1zd79mw+//xzSktLr2p/NDU1YTAYsFgseHl5XbQm6OltV71/u1Ptyil22a744VzNv78umo48amtrrd7Zu7S2tvL1119f07qSkpKYMmUKEydOtBqvqamhoaGB6OhodUyv1xMVFcXu3bsBqKiooL293arGaDQSFham1pSWlmIwGNTgALj77rsxGAxWNWFhYWpwAMTExNDa2kpFRYVaExUVpQZHV83x48epra1Va87vpatm7969tLe3X/T1t7a20tTUZPUnxI3A+VqK8/Ly1P/9wQcfYDAY1McdHR0UFxcTFBR01evLycnhs88+o7y83GZZQ0MDAP7+/lbj/v7+fPXVV2qNq6sr3t7eNjVdz29oaMDPz89m/X5+flY1F27H29sbV1dXq5oLX1vXcxoaGggODr7oevz9/Tl79izffPMN/fv3t+kjMzOT5cuX24wL4eiuKTymTZsGgE6nY8aMGVbLXFxcCAoK4oUXXriqddXV1TFv3jwKCwtxc3O7ZJ1Op7N6rCiKzdiFLqy5WH131HR94rvWmvOlp6eTlpamPm5qaiIwMPCitUI4kmsKj87OTgCCg4MpLy/H19dX84YrKipobGwkPDxcHevo6GDnzp1kZWVx6NAh4Ny7+vnv2I2Njeq7e0BAAG1tbZjNZqujj8bGRkaNGqXWnDhxwmb7J0+etFrPp59+arXcbDbT3t5uVdN1FHL+doAr1jg7O9ucH+qi1+utPgoJcaPQdM6jpqbmuoIDYMKECezfv5+qqir1LyIigkcffZSqqioGDRpEQEAARUVF6nPa2tooKSlRgyE8PBwXFxermvr6ekwmk1oTGRmJxWJhz549as2nn36KxWKxqjGZTNTX16s1hYWF6PV6NdwiIyPZuXOn1de3hYWFGI1G9eNMZGSkVS9dNREREbi4uFzX/hLC0VzTkcf5iouLKS4uprGxUT0i6fL6669f8fmenp6EhYVZjXl4eODj46OOp6amkpGRweDBgxk8eDAZGRn07t2buLg4AAwGAzNnzmTBggX4+PjQt29fFi5cyNChQ9UTsEOGDGHy5MkkJiby8ssvA/C73/2O2NhYQkJCAIiOjiY0NJT4+Hj+/Oc/8+2337Jw4UISExPVM85xcXEsX76chIQEFi9ezBdffEFGRgZ/+MMf1I8ks2fPJisri7S0NBITEyktLeW1117jnXfe0bqbhXBYmsJj+fLlrFixgoiICPr373/FcxBaPfXUU7S0tDBnzhzMZjMjR46ksLAQT09PtWbNmjU4OzvzyCOP0NLSwoQJE9i8eTNOTk5qzdtvv01KSor6TcjUqVPJyspSlzs5ObFt2zbmzJnD6NGjcXd3Jy4ujlWrVqk1BoOBoqIikpKSiIiIwNvbm7S0NKvzFcHBweTn5zN//nw2bNiA0Whk3bp1PPTQQz2yf4SwJ03XefTv35/nn3/e6kIr0T3kOg9hTz1+nUdbW5t6vkAIcXPSFB5PPPEEW7Zs6e5ehBA3EE3nPL7//ns2btzIjh07GDZsmM03Ceff7yGE+HHSFB779u3jzjvvBMBkMlkt66mTp0IIx6IpPD766KPu7kMIcYNxiFvyhRA3Hk1HHuPGjbvsx5MPP/xQc0NCiBuDpvDoOt/Rpb29naqqKkwmk80Nc0KIHydN4bFmzZqLji9btozTp09fV0NCiBtDt57zeOyxx67qvhYhxI2vW8OjtLT0snNzCCF+PDR9bHnwwQetHiuKQn19PXv37mXJkiXd0pgQwrFpCo/zpx8E6NWrFyEhIaxYscJmDk8hxI+TpvDYtGlTd/chhLjBaJ4MCM5NJVhdXY1OpyM0NJThw4d3V19CCAenKTwaGxv59a9/zccff0yfPn1QFAWLxcK4cePIycmhX79+3d2nEMLBaPq2JTk5maamJg4cOMC3336L2WzGZDLR1NRESkpKd/cohHBAmo48CgoK2LFjB0OGDFHHQkND2bBhg5wwFeImoenIo7Oz86Kzgbu4uNhMhiyE+HHSFB7jx49n3rx5HD9+XB37+uuvmT9/PhMmTOi25oQQjktTeGRlZdHc3ExQUBC33XYbP/3pTwkODqa5uZn169d3d49CCAek6ZxHYGAgn332GUVFRRw8eBBFUQgNDbX5sWohxI/XNR15fPjhh4SGhqq/5D5p0iSSk5NJSUlhxIgR3HHHHXzyySc90qgQwrFcU3isXbvW6lfUzmcwGJg1a5ZMfizETeKawuPzzz9n8uTJl1weHR1NRUXFdTclhHB81xQeJ06cuOwPNjs7O3Py5MnrbkoI4fiuKTwGDBjA/v37L7l837599O/f/7qbEkI4vmsKj/vuu48//OEPfP/99zbLWlpaWLp0KbGxsd3WnBDCcV3TV7X/93//x3vvvcfPfvYz5s6dS0hICDqdjurqajZs2EBHRwfPPPNMT/UqhHAg1xQe/v7+7N69myeffJL09HQURQHO/UpcTEwML774Iv7+/j3SqBDCsVzzRWIDBw4kPz8fs9nMkSNHUBSFwYMH4+3t3RP9CSEclObJgLy9vRkxYkR39iKEuIHIz00KITSR8BBCaGLX8MjMzGTEiBF4enri5+fHtGnTOHTokFWNoigsW7YMo9GIu7s7Y8eO5cCBA1Y1ra2tJCcn4+vri4eHB1OnTuXYsWNWNWazmfj4eAwGAwaDgfj4eE6dOmVVc/ToUe6//348PDzw9fUlJSWFtrY2q5r9+/cTFRWFu7s7AwYMYMWKFeqJ4y4lJSWEh4fj5ubGoEGDeOmll65zTwnheOwaHiUlJSQlJVFWVkZRURFnz54lOjqaM2fOqDXPP/88q1evJisri/LycgICApg0aRLNzc1qTWpqKrm5ueTk5LBr1y5Onz5NbGwsHR0dak1cXBxVVVUUFBRQUFBAVVUV8fHx6vKOjg6mTJnCmTNn2LVrFzk5OWzdupUFCxaoNU1NTUyaNAmj0Uh5eTnr169n1apVVvfz1NTUcN999zFmzBgqKytZvHgxKSkpbN26tad2oxB2oVMufNu0o5MnT+Ln50dJSQm/+MUvUBQFo9FIamoqixYtAs4dZfj7+/Pcc88xa9YsLBYL/fr1480332T69OkAHD9+nMDAQPLz84mJiaG6uprQ0FDKysoYOXIkAGVlZURGRnLw4EFCQkLYvn07sbGx1NXVYTQaAcjJySEhIYHGxka8vLzIzs4mPT2dEydOoNfrAVi5ciXr16/n2LFj6HQ6Fi1aRF5eHtXV1errmj17Np9//jmlpaVX3AdNTU0YDAYsFstFb0AECHp6m/adfB1qV06xy3bFD+dq/v11cahzHhaLBYC+ffsC597FGxoarOZF1ev1REVFsXv3buDczz+0t7db1RiNRsLCwtSa0tJSDAaDGhwAd999NwaDwaomLCxMDQ6AmJgYWltb1Zv9SktLiYqKUoOjq+b48ePU1taqNRfO4xoTE8PevXtpb2+/vh0khANxmPBQFIW0tDTuuecewsLCAGhoaACwufDM399fXdbQ0ICrq6vNdSYX1vj5+dls08/Pz6rmwu14e3vj6up62Zqux1eqOXv2LN98841ND62trTQ1NVn9CXEjcJjwmDt3Lvv27eOdd96xWabT6aweK4piM3ahC2suVt8dNedfZXstNV0yMzPVk7gGg4HAwMDLvi4hHIVDhEdycjJ5eXl89NFH3Hrrrep4QEAA8L939S6NjY3qu3tAQABtbW2YzebL1pw4ccJmuydPnrSquXA7ZrOZ9vb2y9Y0NjYCXLHG2dkZHx8fmx7S09OxWCzqX11dnU2NEI7IruGhKApz587lvffe48MPPyQ4ONhqeXBwMAEBARQVFaljbW1tlJSUMGrUKADCw8NxcXGxqqmvr8dkMqk1kZGRWCwW9uzZo9Z8+umnWCwWqxqTyUR9fb1aU1hYiF6vJzw8XK3ZuXOn1de3hYWFGI1GgoKC1Jrze+mqiYiIuOhcKHq9Hi8vL6s/IW4Edg2PpKQk3nrrLbZs2YKnpycNDQ00NDTQ0tICnDvMT01NJSMjg9zcXEwmEwkJCfTu3Zu4uDjg3PSHM2fOZMGCBRQXF1NZWcljjz3G0KFD1QmZhwwZwuTJk0lMTKSsrIyysjISExOJjY0lJCQEODcLWmhoKPHx8VRWVlJcXMzChQutpl2Mi4tDr9eTkJCAyWQiNzeXjIwM0tLS1I8ks2fP5quvviItLY3q6mpef/11XnvtNRYuXPhD714hetR1/dD19crOzgZg7NixVuObNm0iISEBgKeeeoqWlhbmzJmD2Wxm5MiRFBYW4unpqdavWbMGZ2dnHnnkEVpaWpgwYQKbN2/GyclJrXn77bdJSUlRvwmZOnUqWVlZ6nInJye2bdvGnDlzGD16NO7u7sTFxbFq1Sq1xmAwUFRURFJSEhEREXh7e5OWlkZaWppaExwcTH5+PvPnz2fDhg0YjUbWrVvHQw891G37TQhH4FDXeQi5zkPY1w17nYcQ4sYh4SGE0ETCQwihiYSHEEITCQ8hhCYSHkIITSQ8hBCaSHgIITSR8BBCaCLhIYTQRMJDCKGJhIcQQhMJDyGEJhIeQghNJDyEEJpIeAghNJHwEEJoIuEhhNBEwkMIoYmEhxBCEwkPIYQmEh5CCE0kPIQQmkh4CCE0kfAQQmgi4SGE0ETCQwihiYSHEEITCQ8hhCYSHkIITSQ8hBCaSHgIITRxtncD4scj6Oltdtt27copdtv2zUqOPIQQmkh4CCE0kfDoIS+++CLBwcG4ubkRHh7OJ598Yu+WhOhWEh494N133yU1NZVnnnmGyspKxowZw7333svRo0ft3ZoQ3UbCowesXr2amTNn8sQTTzBkyBDWrl1LYGAg2dnZ9m5NiG4j37Z0s7a2NioqKnj66aetxqOjo9m9e7dNfWtrK62trepji8UCQFNT0yW30dn6XTd1e20u1xPYry+4fG9hSz/4ATuxZloeY7dta9G1HxVFuWKthEc3++abb+jo6MDf399q3N/fn4aGBpv6zMxMli9fbjMeGBjYYz1qZVhr7w4uzVF7c9S+rqS5uRmDwXDZGgmPHqLT6aweK4piMwaQnp5OWlqa+rizs5Nvv/0WHx+fi9Zfj6amJgIDA6mrq8PLy6tb1309HLUvcNzeeqovRVFobm7GaDResVbCo5v5+vri5ORkc5TR2NhoczQCoNfr0ev1VmN9+vTpyRbx8vJyqP8QujhqX+C4vfVEX1c64ugiJ0y7maurK+Hh4RQVFVmNFxUVMWrUKDt1JUT3kyOPHpCWlkZ8fDwRERFERkayceNGjh49yuzZs+3dmhDdRsKjB0yfPp3//ve/rFixgvr6esLCwsjPz2fgwIF27Uuv17N06VKbj0n25qh9geP25gh96ZSr+U5GCCEuIOc8hBCaSHgIITSR8BBCaCLhIYTQRMLjJrBz507uv/9+jEYjOp2O999/394tAecuzR8xYgSenp74+fkxbdo0Dh06ZO+2yM7OZtiwYeoFWJGRkWzfvt3ebdnIzMxEp9ORmppql+1LeNwEzpw5w89//nOysrLs3YqVkpISkpKSKCsro6ioiLNnzxIdHc2ZM2fs2tett97KypUr2bt3L3v37mX8+PE88MADHDhwwK59na+8vJyNGzcybNgw+zWhiJsKoOTm5tq7jYtqbGxUAKWkpMTerdjw9vZWXn31VXu3oSiKojQ3NyuDBw9WioqKlKioKGXevHl26UOOPITD6JqOoG/fvnbu5H86OjrIycnhzJkzREZG2rsdAJKSkpgyZQoTJ060ax9yhalwCIqikJaWxj333ENYWJi922H//v1ERkby/fffc8stt5Cbm0toaKi92yInJ4fPPvuM8vJye7ci4SEcw9y5c9m3bx+7du2ydysAhISEUFVVxalTp9i6dSszZsygpKTErgFSV1fHvHnzKCwsxM3NzW59dJHL028yOp2O3Nxcpk2bZu9WVMnJybz//vvs3LmT4OBge7dzURMnTuS2227j5ZdftlsP77//Pr/85S9xcnJSxzo6OtDpdPTq1YvW1larZT1NjjyE3SiKQnJyMrm5uXz88ccOGxxwrtfzp4u0hwkTJrB//36rsd/+9rfcfvvtLFq06AcNDpDwuCmcPn2aI0eOqI9ramqoqqqib9++/OQnP7FbX0lJSWzZsoV//vOfeHp6qhMoGQwG3N3d7dbX4sWLuffeewkMDKS5uZmcnBw+/vhjCgoK7NYTgKenp835IA8PD3x8fOxznsgu3/GIH9RHH32kADZ/M2bMsGtfF+sJUDZt2mTXvh5//HFl4MCBiqurq9KvXz9lwoQJSmFhoV17uhR7flUr5zyEEJrIdR5CCE0kPIQQmkh4CCE0kfAQQmgi4SGE0ETCQwihiYSHEEITCQ8hrsPHH3+MTqfj1KlT9m7lByfhcRNLSEhAp9Oh0+lwdnbmJz/5CU8++SRms/mq11FbW4tOp6OqqqrnGhUOScLjJjd58mTq6+upra3l1Vdf5V//+hdz5syxSy9tbW122a4juhH2hYTHTU6v1xMQEMCtt95KdHQ006dPp7Cw0Kpm06ZNDBkyBDc3N26//XZefPFFdVnXnbDDhw9Hp9MxduxYAMaOHWszMe+0adNISEhQHwcFBfHHP/6RhIQEDAYDiYmJbN68mT59+vDBBx8wZMgQbrnlFjXgLqXro0NxcTERERH07t2bUaNGWU2mnJCQYDMNQWpqqtpvV8/Jycmkpqbi7e2Nv78/Gzdu5MyZM/z2t7/F09OT22677aKTIf/73//m5z//OW5ubowcOdLm7tfdu3fzi1/8And3dwIDA0lJSbGaq/Vi+6KtrY25c+fSv39/3NzcCAoKIjMz85L74Ycm4SFUX375JQUFBbi4uKhjr7zyCs888wx/+tOfqK6uJiMjgyVLlvDXv/4VgD179gCwY8cO6uvree+9965pm3/+858JCwujoqKCJUuWAPDdd9+xatUq3nzzTXbu3MnRo0dZuHDhFdf1zDPP8MILL7B3716cnZ15/PHHr6kXgL/+9a/4+vqyZ88ekpOTefLJJ3n44YcZNWoUn332GTExMcTHx/Pdd99ZPe/3v/89q1atory8HD8/P6ZOnUp7eztwblaymJgYHnzwQfbt28e7777Lrl27mDt37mX3xbp168jLy+Nvf/sbhw4d4q233iIoKOiaX1OPscvteMIhzJgxQ3FyclI8PDwUNzc39a7W1atXqzWBgYHKli1brJ737LPPKpGRkYqiKEpNTY0CKJWVlVY1F7vb84EHHrC6k3fgwIHKtGnTrGo2bdqkAMqRI0fUsQ0bNij+/v6XfB1ddw3v2LFDHdu2bZsCKC0tLeprfeCBB6yeN2/ePCUqKsqq53vuuUd9fPbsWcXDw0OJj49Xx+rr6xVAKS0ttdp2Tk6OWvPf//5XcXd3V959911FURQlPj5e+d3vfme17U8++UTp1auX2t/F9kVycrIyfvx4pbOz85Kv3Z5kPo+b3Lhx48jOzua7777j1Vdf5fDhwyQnJwNw8uRJ6urqmDlzJomJiepzzp49i8Fg6JbtR0RE2Iz17t2b2267TX3cv39/Ghsbr7iu83+GoH///gA0NjZe05wl56/DyckJHx8fhg4dqo75+/ur6z3f+ZMj9+3bl5CQEKqrqwGoqKjgyJEjvP3222qNoih0dnZSU1PDkCFDANt9kZCQwKRJkwgJCWHy5MnExsYSHR191a+lp0l43OQ8PDz46U9/CsC6desYN24cy5cv59lnn6WzsxM499Fl5MiRVs+70qxVvXr1Qrlgtoeuw/gLt3+h8z82wbmpEy9c18Wc/zydTgegvoar7edi277cei/n/NpZs2aRkpJiU3N+sF24L+666y5qamrYvn07O3bs4JFHHmHixIn84x//uOK2fwgSHsLK0qVLuffee3nyyScxGo0MGDCAL7/8kkcfffSi9a6ursC5uTTP169fP6uTnB0dHZhMJsaNG9dzzV9Gv379MJlMVmNVVVU2YaFVWVmZGgRms5nDhw9z++23A+dC4MCBA2pIXwsvLy+mT5/O9OnT+dWvfsXkyZP59ttvHeLnKeSEqbAyduxY7rjjDjIyMgBYtmwZmZmZ/OUvf+Hw4cPs37+fTZs2sXr1agD8/Pxwd3enoKCAEydOqL+9Mn78eLZt28a2bds4ePAgc+bMseuFVOPHj2fv3r288cYbfPHFFyxdutQmTK7HihUrKC4uxmQykZCQgK+vr/rtzqJFiygtLSUpKYmqqiq++OIL8vLy1I+Hl7JmzRpycnI4ePAghw8f5u9//zsBAQH06dOn2/q+HhIewkZaWhqvvPIKdXV1PPHEE7z66qts3ryZoUOHEhUVxebNm9WvaJ2dnVm3bh0vv/wyRqORBx54AIDHH3+cGTNm8Jvf/IaoqCiCg4PtdtQBEBMTw5IlS3jqqacYMWIEzc3N/OY3v+m29a9cuZJ58+YRHh5OfX09eXl56lHZsGHDKCkp4YsvvmDMmDEMHz6cJUuWqOdlLuWWW27hueeeIyIighEjRlBbW0t+fj69ejnGf7YyDaEQQhPHiDAhxA1HwkMIoYmEhxBCEwkPIYQmEh5CCE0kPIQQmkh4CCE0kfAQQmgi4SGE0ETCQwihiYSHEEITCQ8hhCb/DznOO5/JgFxPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 800x200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 4656322\n",
      "2 225862\n",
      "3 22117\n",
      "4 396\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQ8AAADVCAYAAABJ/XXDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAnKUlEQVR4nO3dfVRU953H8ffIw4AERhQBR0nUtIdI0TQBazCx+ASYiA/b7OoeIpXGsKIIItpQzHFjbAtqqaY+xyZq18Sd7K4lpUekIIkYVkYBQRkX83CKAcuTJiMooYB49w8Pt47gAzeYwfh9nTPndO79zr3fe2s+3IeZ+9MpiqIghBC9NMDeDQghHkwSHkIITSQ8hBCaSHgIITSR8BBCaCLhIYTQRMJDCKGJhIcQQhNHezcgbF2/fp3a2lrc3d3R6XT2bkc8ZBRF4cqVKxiNRgYMuPOxhYRHP1NbW4ufn5+92xAPuZqaGkaMGHHHGgmPfsbd3R248X+eh4eHnbsRD5vm5mb8/PzUf4d3IuHRz3Sdqnh4eEh4CLu5l1NmuWAqhNBEwkMIoYmEhxBCEwkPIYQmcsH0ATTyF4fsst7z62faZb2if5IjDyGEJhIeQghNJDyEEJpIeAghNJHwEEJoIuEhhNBEwkMIoYmEhxBCEwkPIYQmEh5CCE0kPIQQmkh4CCE0kfAQQmgi4SGE0ETCQwihSb8Jj/T0dHQ6HUlJSeo0RVFYu3YtRqMRV1dXJk+ezNmzZ20+19bWRkJCAl5eXri5uTF79mwuXLhgU2O1WomOjsZgMGAwGIiOjuby5cs2NdXV1cyaNQs3Nze8vLxITEykvb3dpqaiooLQ0FBcXV0ZPnw469atQ1EUm5qCggKCgoJwcXFh9OjR7Nq165vvHCH6oX4RHsXFxezevZtx48bZTN+4cSObNm1i27ZtFBcX4+vrS1hYGFeuXFFrkpKSyMzMxGQyUVhYyNWrV4mMjKSzs1OtiYqKory8nJycHHJycigvLyc6Olqd39nZycyZM2lpaaGwsBCTycTBgwdZuXKlWtPc3ExYWBhGo5Hi4mK2bt1KRkYGmzZtUmuqqqp44YUXmDRpEmVlZaxevZrExEQOHjx4P3abEHalU2790/ktu3r1Kk8//TQ7duzgV7/6FT/84Q958803URQFo9FIUlISKSkpwI2jDB8fHzZs2MDixYtpampi6NCh7N+/n/nz5wP/GDQpOzubiIgIKisrCQgIwGw2M2HCBADMZjMhISGcO3cOf39/Dh8+TGRkJDU1NRiNRgBMJhMxMTE0Njbi4eHBzp07SU1NpaGhAb1eD8D69evZunUrFy5cQKfTkZKSQlZWFpWVler2xcXFcfr0aYqKiu5pfzQ3N2MwGGhqarrt0AvyJDFxv9zLv78udj/yiI+PZ+bMmUyfPt1melVVFfX19YSHh6vT9Ho9oaGhHD9+HIDS0lI6OjpsaoxGI4GBgWpNUVERBoNBDQ6AZ555BoPBYFMTGBioBgdAREQEbW1tlJaWqjWhoaFqcHTV1NbWcv78ebXm5l66akpKSujo6Ohx+9va2mhubrZ5CfEgsGt4mEwmTp06RXp6erd59fX1APj4+NhM9/HxUefV19fj7OyMp6fnHWu8vb27Ld/b29um5tb1eHp64uzsfMearvd3q7l27RqXLl3qcR+kp6er12IMBoMMNSkeGHYLj5qaGpYvX867776Li4vLbetuHblKUZS7jmZ1a01P9X1R03XG19uam6WmptLU1KS+ampqbrtdQvQndguP0tJSGhsbCQoKwtHREUdHRwoKCtiyZQuOjo7d/qp3aWxsVOf5+vrS3t6O1Wq9Y01DQ0O39V+8eNGm5tb1WK1WOjo67ljT2NgIcNcaR0dHhgwZ0uN+0Ov16tCSMsSkeJDYLTymTZtGRUUF5eXl6is4OJiXXnqJ8vJyRo8eja+vL3l5eepn2tvbKSgoYOLEiQAEBQXh5ORkU1NXV4fFYlFrQkJCaGpq4uTJk2rNiRMnaGpqsqmxWCzU1dWpNbm5uej1eoKCgtSaY8eO2dy+zc3NxWg0MnLkSLXm5l66aoKDg3FycuqL3SZEv2G38HB3dycwMNDm5ebmxpAhQwgMDFS/85GWlkZmZiYWi4WYmBgGDhxIVFQUAAaDgUWLFrFy5Ury8/MpKytjwYIFjB07Vr0AO2bMGGbMmEFsbCxmsxmz2UxsbCyRkZH4+/sDEB4eTkBAANHR0ZSVlZGfn8+qVauIjY1VjwSioqLQ6/XExMRgsVjIzMwkLS2N5ORk9ZQkLi6OL774guTkZCorK9mzZw/vvPMOq1atssMeFuL+6teDPr366qu0traydOlSrFYrEyZMIDc3F3d3d7Vm8+bNODo6Mm/ePFpbW5k2bRr79u3DwcFBrXnvvfdITExU74TMnj2bbdu2qfMdHBw4dOgQS5cu5dlnn8XV1ZWoqCgyMjLUGoPBQF5eHvHx8QQHB+Pp6UlycjLJyclqzahRo8jOzmbFihVs374do9HIli1bePHFF+/nbhLCLuz+PQ9hS77nIezpgfqehxDiwSThIYTQRMJDCKGJhIcQQhMJDyGEJhIeQghNJDyEEJpIeAghNJHwEEJoIuEhhNBEwkMIoYmEhxBCEwkPIYQmEh5CCE0kPIQQmkh4CCE0kfAQQmgi4SGE0ETCQwihiV3DY+fOnYwbN04dryQkJITDhw+r8xVFYe3atRiNRlxdXZk8eTJnz561WUZbWxsJCQl4eXnh5ubG7NmzuXDhgk2N1WolOjpaHZUtOjqay5cv29RUV1cza9Ys3Nzc8PLyIjEx0WaYBYCKigpCQ0NxdXVl+PDhrFu3jlsfAVtQUEBQUBAuLi6MHj2aXbt29cGeEqL/sWt4jBgxgvXr11NSUkJJSQlTp05lzpw5akBs3LiRTZs2sW3bNoqLi/H19SUsLIwrV66oy0hKSiIzMxOTyURhYSFXr14lMjKSzs5OtSYqKory8nJycnLIycmhvLyc6OhodX5nZyczZ86kpaWFwsJCTCYTBw8eZOXKlWpNc3MzYWFhGI1GiouL2bp1KxkZGWzatEmtqaqq4oUXXmDSpEmUlZWxevVqEhMTOXjw4P3cjULYRb97evrgwYP5zW9+w8svv4zRaCQpKYmUlBTgxlGGj48PGzZsYPHixTQ1NTF06FD279/P/PnzAaitrcXPz4/s7GwiIiKorKwkICAAs9msDnZtNpsJCQnh3Llz+Pv7c/jwYSIjI6mpqVEHuzaZTMTExNDY2IiHhwc7d+4kNTWVhoYGdbDr9evXs3XrVi5cuIBOpyMlJYWsrCwqKyvV7YmLi+P06dMUFRXd0/bL09OFPT2QT0/v7OzEZDLR0tJCSEgIVVVV1NfX24w6r9frCQ0NVUe3Ly0tpaOjw6bGaDQSGBio1hQVFWEwGNTgAHjmmWcwGAw2NYGBgWpwwI3R7dva2igtLVVrQkND1eDoqqmtreX8+fNqzc29dNWUlJTQ0dHRF7tJiH5DU3iMHj2aL7/8stv0y5cvM3r06F4tq6KigkceeQS9Xk9cXByZmZkEBASoY772NOr8zaPSOzs74+npeccab2/vbuv19va+4+j2np6eODs737Hm1vF0b1dz7do1Ll261OP2t7W10dzcbPMS4kGgKTzOnz9vc02hS1tbG3/72996tSx/f3/Ky8sxm80sWbKEhQsX8n//93/q/J5Gnb/diPO3q+mpvi9qus74eltzs/T0dPVCrsFgwM/P77bbJUR/0qvhJrOystT//Ze//AWDwaC+7+zsJD8/Xx30+V45Ozvzve99D4Dg4GCKi4v53e9+p17nqK+vZ9iwYWp9Y2Ojzaj07e3tWK1Wm6OPxsZGdRBrX19fGhoauq334sWLNss5ceKEzXyr1UpHR4dNTdcRxs3rAe5a4+joyJAhQ3rc/tTUVJshK5ubmyVAxAOhV+Exd+5c4MZf0YULF9rMc3JyYuTIkfz2t7/9Rg0pikJbWxujRo3C19eXvLw8nnrqKQDa29spKChgw4YNAAQFBeHk5EReXh7z5s0DoK6uDovFwsaNG4EbI9c3NTVx8uRJfvSjHwFw4sQJmpqa1IAJCQnh17/+NXV1dWpQ5ebmotfrCQoKUmtWr15Ne3s7zs7Oao3RaFQDMyQkhD//+c8225Obm0twcDBOTk49bq9er7e5jiLEg6JXpy3Xr1/n+vXrPProozQ2Nqrvr1+/TltbG5988gmRkZH3vLzVq1fz8ccfc/78eSoqKnjttdc4evQoL730EjqdjqSkJNLS0sjMzMRisRATE8PAgQOJiooCbgw+vWjRIlauXEl+fj5lZWUsWLCAsWPHMn36dADGjBnDjBkziI2NxWw2YzabiY2NJTIyEn9/fwDCw8MJCAggOjqasrIy8vPzWbVqFbGxseoV56ioKPR6PTExMVgsFjIzM0lLSyM5OVk9JYmLi+OLL74gOTmZyspK9uzZwzvvvMOqVat6s5uFeCD06sijS1VVVZ+svKGhgejoaOrq6jAYDIwbN46cnBzCwsIAePXVV2ltbWXp0qVYrVYmTJhAbm4u7u7u6jI2b96Mo6Mj8+bNo7W1lWnTprFv3z4cHBzUmvfee4/ExET1Tsjs2bPZtm2bOt/BwYFDhw6xdOlSnn32WVxdXYmKiiIjI0OtMRgM5OXlER8fT3BwMJ6eniQnJ9uccowaNYrs7GxWrFjB9u3bMRqNbNmyhRdffLFP9pcQ/Ynm73nk5+eTn5+vHoHcbM+ePX3S3MNIvuch7Kk33/PQdOTxxhtvsG7dOoKDgxk2bNhd734IIb57NIXHrl272Ldvn81XvIUQDxdN3/Nob29X71QIIR5OmsLjlVde4cCBA33dixDiAaLptOXvf/87u3fv5siRI4wbN67bdxhu/qWpEOK7SVN4nDlzhh/+8IcAWCwWm3ly8VSIh4Om8Pjoo4/6ug8hxAOm3/wkXwjxYNF05DFlypQ7np58+OGHmhsSQjwYNIVH1/WOLh0dHZSXl2OxWLr9YE4I8d2kKTw2b97c4/S1a9dy9erVb9SQEOLB0KfXPBYsWCC/axHiIdGn4VFUVISLi0tfLlII0U9pOm35yU9+YvNeURTq6uooKSlhzZo1fdKYEKJ/0xQeNz9+EGDAgAH4+/uzbt26bk8PF0J8N2kKj7179/Z1H0KIB4ym8OhSWlpKZWUlOp2OgIAA9VmjQojvPk3h0djYyL/+679y9OhRBg0ahKIoNDU1MWXKFEwmE0OHDu3rPoUQ/Yymuy0JCQk0Nzdz9uxZvvrqK6xWKxaLhebmZhITE/u6RyFEP6TpyCMnJ4cjR44wZswYdVpAQADbt2+XC6ZCPCQ0HXlcv369x3FInJycuj0M+U7S09MZP3487u7ueHt7M3fuXD755BObGkVRWLt2LUajEVdXVyZPnszZs2dtatra2khISMDLyws3Nzdmz57NhQsXbGqsVivR0dHqyGzR0dFcvnzZpqa6uppZs2bh5uaGl5cXiYmJtLe329RUVFQQGhqKq6srw4cPZ926ddz6DOmCggKCgoJwcXFh9OjR7Nq16573iRAPCk3hMXXqVJYvX05tba067W9/+xsrVqxg2rRp97ycgoIC4uPjMZvN5OXlce3aNcLDw2lpaVFrNm7cyKZNm9i2bRvFxcX4+voSFhbGlStX1JqkpCQyMzMxmUwUFhZy9epVIiMjbYbEjIqKory8nJycHHJycigvL7d5BmtnZyczZ86kpaWFwsJCTCYTBw8eZOXKlWpNc3MzYWFhGI1GiouL2bp1KxkZGTYPP6qqquKFF15g0qRJlJWVsXr1ahITEzl48OC972AhHgCahl6oqalhzpw5WCwW/Pz80Ol0VFdXM3bsWP70pz8xYsQITc1cvHgRb29vCgoK+PGPf4yiKBiNRpKSktThJ9va2vDx8WHDhg0sXryYpqYmhg4dyv79+5k/fz4AtbW1+Pn5kZ2dTUREBJWVlQQEBGA2m5kwYQIAZrOZkJAQzp07h7+/P4cPHyYyMpKamhqMRiMAJpOJmJgYGhsb8fDwYOfOnaSmptLQ0KCO8rZ+/Xq2bt3KhQsX0Ol0pKSkkJWVRWVlpbpdcXFxnD59mqKiorvuAxl6QdhTb4Ze0HTk4efnx6lTpzh06BBJSUkkJiaSnZ1NaWmp5uAAaGpqAmDw4MHAjb/i9fX1NtdR9Ho9oaGhHD9+HLhxu7ijo8Omxmg0EhgYqNYUFRVhMBjU4AB45plnMBgMNjWBgYFqcABERETQ1tZGaWmpWhMaGmozPGRERAS1tbWcP39erbn1uk9ERAQlJSV0dHR02+a2tjaam5ttXkI8CHoVHh9++CEBAQHqP/CwsDASEhJITExk/Pjx/OAHP+Djjz/W1IiiKCQnJ/Pcc88RGBgIoA4a3TWQdBcfHx91Xn19Pc7OzjYDXfdU4+3t3W2d3t7eNjW3rsfT0xNnZ+c71nS9v1vNtWvXuHTpUrce0tPT1eswBoNBBrkWD4xehcebb75pM37rzQwGA4sXL9b88ONly5Zx5swZ/vM//7PbvFsfPKQoyl2flXprTU/1fVHTddbX25ouqampNDU1qa+ampo7bpcQ/UWvwuP06dPMmDHjtvPDw8PVQ/zeSEhIICsri48++sjmtMfX1xf4x1/1Lo2Njepfd19fX9rb27FarXesaWho6Lbeixcv2tTcuh6r1UpHR8cdaxobGwHuWuPo6MiQIUO69aDX6/Hw8LB5CfEg6FV4NDQ09HiLtoujoyMXL1685+UpisKyZcv44x//yIcffsioUaNs5o8aNQpfX1/y8vLUae3t7RQUFKiDTgUFBeHk5GRTU1dXh8ViUWtCQkJoamri5MmTas2JEydoamqyqbFYLNTV1ak1ubm56PV6goKC1Jpjx47Z3L7Nzc3FaDQycuRItebmXrpqgoOD77jvhHjQ9Co8hg8fTkVFxW3nnzlzhmHDht3z8uLj43n33Xc5cOAA7u7u1NfXU19fT2trK3DjMD8pKYm0tDQyMzOxWCzExMQwcOBAoqKigBunS4sWLWLlypXk5+dTVlbGggULGDt2LNOnTwdgzJgxzJgxg9jYWMxmM2azmdjYWCIjI/H39wduHDUFBAQQHR1NWVkZ+fn5rFq1yuY0LSoqCr1eT0xMDBaLhczMTNLS0khOTlZPSeLi4vjiiy9ITk6msrKSPXv28M4777Bq1are7Goh+r1e3apNSEjg6NGjFBcXd3voT2trKz/60Y+YMmUKW7ZsubeV3+a6xd69e4mJiQFuHJ288cYbvPXWW1itViZMmMD27dvVi6pwYxCqn//85xw4cIDW1lamTZvGjh07bC4+fvXVVyQmJpKVlQXA7Nmz2bZtG4MGDVJrqqurWbp0KR9++CGurq5ERUWRkZFhc3eloqKC+Ph4Tp48iaenJ3Fxcfz7v/+7zbYUFBSwYsUKzp49i9FoJCUlhbi4uHvaJ3KrVthTb27V9io8GhoaePrpp3FwcGDZsmX4+/uj0+morKxk+/btdHZ2curUqW53G8S9k/AQ9tSb8OjVb1t8fHw4fvw4S5YsITU11eYuQkREBDt27JDgEOIh0esfxj322GNkZ2djtVr5/PPPURSF73//+92+ZyGE+G7T/DAgT09Pxo8f35e9CCEeIDLcpBBCEwkPIYQmEh5CCE0kPIQQmkh4CCE0kfAQQmgi4SGE0ETCQwihiYSHEEITCQ8hhCYSHkIITSQ8hBCaSHgIITSR8BBCaCLhIYTQRMJDCKGJhIcQQhO7hsexY8eYNWsWRqMRnU7HBx98YDNfURTWrl2L0WjE1dWVyZMnc/bsWZuatrY2EhIS8PLyws3NjdmzZ3PhwgWbGqvVSnR0tDqkY3R0NJcvX7apqa6uZtasWbi5ueHl5UViYqLN+Cxw48npoaGhuLq6Mnz4cNatW8etz48uKCggKCgIFxcXRo8eza5du77ZThKin7JreLS0tPDkk0+ybdu2Hudv3LiRTZs2sW3bNoqLi/H19SUsLIwrV66oNUlJSWRmZmIymSgsLOTq1atERkbS2dmp1kRFRVFeXk5OTg45OTmUl5cTHR2tzu/s7GTmzJm0tLRQWFiIyWTi4MGDrFy5Uq1pbm4mLCwMo9FIcXExW7duJSMjw2Z4zaqqKl544QUmTZpEWVkZq1evJjExkYMHD/blbhOiX+jV0Av3k06nIzMzk7lz5wI3jjqMRiNJSUmkpKQAN44yfHx82LBhA4sXL6apqYmhQ4eyf/9+5s+fD0BtbS1+fn5kZ2cTERFBZWUlAQEBmM1mJkyYAIDZbCYkJIRz587h7+/P4cOHiYyMpKamBqPRCIDJZCImJobGxkY8PDzYuXMnqampNDQ0qOO4rF+/nq1bt3LhwgV0Oh0pKSlkZWVRWVmpbldcXBynT5+mqKjonvaDDL0g7Kk3Qy/022seVVVV1NfXEx4erk7T6/WEhoZy/PhxAEpLS+no6LCpMRqNBAYGqjVFRUUYDAY1OACeeeYZDAaDTU1gYKAaHAARERG0tbWpY+8WFRURGhpqMwBUREQEtbW1nD9/Xq25uZeumpKSEjo6Onrczra2Npqbm21eQjwI+m14dA0Wfes4MD4+Puq8+vp6nJ2duw37cGuNt7d3t+V7e3vb1Ny6Hk9PT5ydne9Y0/X+bjXXrl3j0qVLPW5nenq6ei3GYDDYjHInRH/Wb8Ojy61DUiqKctthKm9X01N9X9TcPOhVb2pulpqaSlNTk/qqqam57XYJ0Z/02/Dw9fUF/vFXvUtjY6P6193X15f29nasVusdaxoaGrot/+LFizY1t67HarXS0dFxx5rGxkaAu9Y4OjoyZMiQHrdTr9fj4eFh8xLiQdBvw2PUqFH4+vqSl5enTmtvb6egoICJEycCEBQUhJOTk01NXV0dFotFrQkJCaGpqYmTJ0+qNSdOnKCpqcmmxmKxUFdXp9bk5uai1+sJCgpSa44dO2Zz+zY3Nxej0cjIkSPVmpt76aoJDg7GycmpL3aLEP2GXcPj6tWrlJeXU15eDty4SFpeXk51dTU6nY6kpCTS0tLIzMzEYrEQExPDwIEDiYqKAsBgMLBo0SJWrlxJfn4+ZWVlLFiwgLFjxzJ9+nQAxowZw4wZM4iNjcVsNmM2m4mNjSUyMhJ/f38AwsPDCQgIIDo6mrKyMvLz81m1ahWxsbHqkUBUVBR6vZ6YmBgsFguZmZmkpaWRnJysnpLExcXxxRdfkJycTGVlJXv27OGdd95h1apV3/KeFeL+0zzcZF8oKSlhypQp6vvk5GQAFi5cyL59+3j11VdpbW1l6dKlWK1WJkyYQG5uLu7u7upnNm/ejKOjI/PmzaO1tZVp06axb98+HBwc1Jr33nuPxMRE9U7I7Nmzbb5b4uDgwKFDh1i6dCnPPvssrq6uREVFkZGRodYYDAby8vKIj48nODgYT09PkpOT1Z7hxtFSdnY2K1asYPv27RiNRrZs2cKLL77Y9ztPCDvrN9/zEDfI9zyEPX0nvuchhOjfJDyEEJpIeAghNJHwEEJoIuEhhNBEwkMIoYmEhxBCEwkPIYQmEh5CCE0kPIQQmkh4CCE0kfAQQmgi4SGE0ETCQwihiYSHEEITCQ8hhCYSHkIITez6GELx3WKvJ5yBPOXMHuTIQwihiYSHEEITCY/7ZMeOHYwaNQoXFxeCgoL4+OOP7d2SEH1KwuM+eP/990lKSuK1116jrKyMSZMm8fzzz1NdXW3v1oToMxIe98GmTZtYtGgRr7zyCmPGjOHNN9/Ez8+PnTt32rs1IfqM3G3pY+3t7ZSWlvKLX/zCZnp4eDjHjx/vVt/W1kZbW5v6vqmpCbgxfsbtXG/7uo+67Z079QT26wvu3Fvg63/5FjuxZXkjwm7r1qJrP97LcE4SHn3s0qVLdHZ2qoNfd/Hx8ek2CDZAeno6b7zxRrfpfn5+961HrQxv2ruD2+uvvfXXvu7mypUrGAyGO9ZIeNwnXePXdlEUpds0gNTUVJshK69fv85XX33FkCFDeqz/Jpqbm/Hz86Ompuauo4F9m/prX9B/e7tffSmKwpUrVzAajXetlfDoY15eXjg4OHQ7ymhsbOx2NAKg1+vR6/U20wYNGnQ/W8TDw6Nf/YfQpb/2Bf23t/vR192OOLrIBdM+5uzsTFBQEHl5eTbT8/LymDhxop26EqLvyZHHfZCcnEx0dDTBwcGEhISwe/duqquriYuLs3drQvQZCY/7YP78+Xz55ZesW7eOuro6AgMDyc7O5rHHHrNrX3q9ntdff73baZK99de+oP/21h/60in3ck9GCCFuIdc8hBCaSHgIITSR8BBCaCLhIYTQRMLjIXDs2DFmzZqF0WhEp9PxwQcf2Lsl4MZX88ePH4+7uzve3t7MnTuXTz75xN5tsXPnTsaNG6d+ASskJITDhw/bu61u0tPT0el0JCUl2WX9Eh4PgZaWFp588km2bdtm71ZsFBQUEB8fj9lsJi8vj2vXrhEeHk5LS4td+xoxYgTr16+npKSEkpISpk6dypw5czh79qxd+7pZcXExu3fvZty4cfZrQhEPFUDJzMy0dxs9amxsVACloKDA3q104+npqbz99tv2bkNRFEW5cuWK8v3vf1/Jy8tTQkNDleXLl9ulDznyEP1G1+MIBg8ebOdO/qGzsxOTyURLSwshISH2bgeA+Ph4Zs6cyfTp0+3ah3zDVPQLiqKQnJzMc889R2BgoL3boaKigpCQEP7+97/zyCOPkJmZSUBAgL3bwmQycerUKYqLi+3dioSH6B+WLVvGmTNnKCwstHcrAPj7+1NeXs7ly5c5ePAgCxcupKCgwK4BUlNTw/Lly8nNzcXFxcVufXSRr6c/ZHQ6HZmZmcydO9feragSEhL44IMPOHbsGKNGjbJ3Oz2aPn06jz/+OG+99Zbdevjggw/4p3/6JxwcHNRpnZ2d6HQ6BgwYQFtbm828+02OPITdKIpCQkICmZmZHD16tN8GB9zo9ebHRdrDtGnTqKiosJn2s5/9jCeeeIKUlJRvNThAwuOhcPXqVT7//HP1fVVVFeXl5QwePJhHH33Ubn3Fx8dz4MAB/vSnP+Hu7q4+QMlgMODq6mq3vlavXs3zzz+Pn58fV65cwWQycfToUXJycuzWE4C7u3u360Fubm4MGTLEPteJ7HKPR3yrPvroIwXo9lq4cKFd++qpJ0DZu3evXft6+eWXlccee0xxdnZWhg4dqkybNk3Jzc21a0+3Y89btXLNQwihiXzPQwihiYSHEEITCQ8hhCYSHkIITSQ8hBCaSHgIITSR8BBCaCLhIcQ3cPToUXQ6HZcvX7Z3K986CY+HWExMDDqdDp1Oh6OjI48++ihLlizBarXe8zLOnz+PTqejvLz8/jUq+iUJj4fcjBkzqKur4/z587z99tv8+c9/ZunSpXbppb293S7r7Y8ehH0h4fGQ0+v1+Pr6MmLECMLDw5k/fz65ubk2NXv37mXMmDG4uLjwxBNPsGPHDnVe1y9hn3rqKXQ6HZMnTwZg8uTJ3R7MO3fuXGJiYtT3I0eO5Fe/+hUxMTEYDAZiY2PZt28fgwYN4i9/+QtjxozhkUceUQPudrpOHfLz8wkODmbgwIFMnDjR5mHKMTEx3R5DkJSUpPbb1XNCQgJJSUl4enri4+PD7t27aWlp4Wc/+xnu7u48/vjjPT4M+X//93958skncXFxYcKECd1+/Xr8+HF+/OMf4+rqip+fH4mJiTbPau1pX7S3t7Ns2TKGDRuGi4sLI0eOJD09/bb74dsm4SFUf/3rX8nJycHJyUmd9vvf/57XXnuNX//611RWVpKWlsaaNWv4wx/+AMDJkycBOHLkCHV1dfzxj3/s1Tp/85vfEBgYSGlpKWvWrAHg66+/JiMjg/3793Ps2DGqq6tZtWrVXZf12muv8dvf/paSkhIcHR15+eWXe9ULwB/+8Ae8vLw4efIkCQkJLFmyhH/5l39h4sSJnDp1ioiICKKjo/n6669tPvfzn/+cjIwMiouL8fb2Zvbs2XR0dAA3nkoWERHBT37yE86cOcP7779PYWEhy5Ytu+O+2LJlC1lZWfzXf/0Xn3zyCe+++y4jR47s9TbdN3b5OZ7oFxYuXKg4ODgobm5uiouLi/qr1k2bNqk1fn5+yoEDB2w+98tf/lIJCQlRFEVRqqqqFEApKyuzqenp155z5syx+SXvY489psydO9emZu/evQqgfP755+q07du3Kz4+Prfdjq5fDR85ckSddujQIQVQWltb1W2dM2eOzeeWL1+uhIaG2vT83HPPqe+vXbumuLm5KdHR0eq0uro6BVCKiops1m0ymdSaL7/8UnF1dVXef/99RVEUJTo6Wvm3f/s3m3V//PHHyoABA9T+etoXCQkJytSpU5Xr16/fdtvtSZ7n8ZCbMmUKO3fu5Ouvv+btt9/m008/JSEhAYCLFy9SU1PDokWLiI2NVT9z7do1DAZDn6w/ODi427SBAwfy+OOPq++HDRtGY2PjXZd18zAEw4YNA6CxsbFXzyy5eRkODg4MGTKEsWPHqtN8fHzU5d7s5ocjDx48GH9/fyorKwEoLS3l888/57333lNrFEXh+vXrVFVVMWbMGKD7voiJiSEsLAx/f39mzJhBZGQk4eHh97wt95uEx0POzc2N733vewBs2bKFKVOm8MYbb/DLX/6S69evAzdOXSZMmGDzubs9tWrAgAEotzztoesw/tb13+rm0ya48ejEW5fVk5s/p9PpANRtuNd+elr3nZZ7JzfXLl68mMTExG41Nwfbrfvi6aefpqqqisOHD3PkyBHmzZvH9OnT+Z//+Z+7rvvbIOEhbLz++us8//zzLFmyBKPRyPDhw/nrX//KSy+91GO9s7MzcONZmjcbOnSozUXOzs5OLBYLU6ZMuX/N38HQoUOxWCw208rLy7uFhVZms1kNAqvVyqeffsoTTzwB3AiBs2fPqiHdGx4eHsyfP5/58+fzz//8z8yYMYOvvvqqXwxPIRdMhY3Jkyfzgx/8gLS0NADWrl1Leno6v/vd7/j000+pqKhg7969bNq0CQBvb29cXV3JycmhoaFBHXtl6tSpHDp0iEOHDnHu3DmWLl1q1y9STZ06lZKSEv7jP/6Dzz77jNdff71bmHwT69atIz8/H4vFQkxMDF5eXurdnZSUFIqKioiPj6e8vJzPPvuMrKws9fTwdjZv3ozJZOLcuXN8+umn/Pd//ze+vr4MGjSoz/r+JiQ8RDfJycn8/ve/p6amhldeeYW3336bffv2MXbsWEJDQ9m3b596i9bR0ZEtW7bw1ltvYTQamTNnDgAvv/wyCxcu5Kc//SmhoaGMGjXKbkcdABEREaxZs4ZXX32V8ePHc+XKFX7605/22fLXr1/P8uXLCQoKoq6ujqysLPWobNy4cRQUFPDZZ58xadIknnrqKdasWaNel7mdRx55hA0bNhAcHMz48eM5f/482dnZDBjQP/6zlccQCiE06R8RJoR44Eh4CCE0kfAQQmgi4SGE0ETCQwihiYSHEEITCQ8hhCYSHkIITSQ8hBCaSHgIITSR8BBCaCLhIYTQ5P8BHFbRexC3uOEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 800x200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Make a plot of return numbers from lasFile\n",
    "plotReturnNumberCount(lasData1)\n",
    "plotReturnNumberCount(lasData2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56b1032d-e974-46f3-a07d-ac0db5d2d8ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select a 100 m square subset and export to a new laz file\n",
    "lasFile = r'E:\\lidar_for_jack\\katoomba_201804_lid1\\260000_6244000_Katoomba201804-LID1-C3-AHD_z56.las'\n",
    "outFile = r'E:\\lidar_for_jack\\katoomba_201804_lid1\\test.laz'\n",
    "epsg = 28355\n",
    "\n",
    "#min and max X: 260000.000000 262000.000000\n",
    "#min and max Y: 6244000.000000 6246000.000000\n",
    "\n",
    "minX = 261000\n",
    "maxX = 261100\n",
    "minY = 6245000\n",
    "maxY = 6245100\n",
    "\n",
    "las = laspy.read(lasFile)\n",
    "subset = laspy.create(point_format=las.header.point_format, file_version=las.header.version)\n",
    "subset.points = las.points[(las.x >= minX) & (las.x <= maxX) & (las.y >= minY) & (las.y <= maxY)]\n",
    "subset.header.scales = las.header.scales\n",
    "subset.header.offset = las.header.offset\n",
    "subset.header.add_crs(crs.CRS.from_user_input(epsg))\n",
    "subset.write(outFile)\n",
    "\n",
    "# Check the new file\n",
    "print('min and max X: %f %f'%(np.min(subset.points.x), np.max(subset.points.x)))\n",
    "print('min and max Y: %f %f'%(np.min(subset.points.y), np.max(subset.points.y)))\n",
    "print('min and max Z: %f %f'%(np.min(subset.points.z), np.max(subset.points.z)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74184f41-0c31-4bed-8be6-bbfac4578a4a",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}