From ae9d85d545135c28ba3327eba4c82e9c47cd3b24 Mon Sep 17 00:00:00 2001 From: Josh von Nonn <68934869+vonnonn@users.noreply.github.com> Date: Fri, 1 Dec 2023 10:52:14 -0800 Subject: [PATCH] Add files via upload This is a notebook that allows the user to iterate ground classifier parameters and visualize results for optimizing the algorithm of their choice using PDAL. Caveat- this notebook requires hard paths and the knowledge to build a Python environment and import packages. This notebook is also designed for SfM RGB point clouds. I will be adding another notebook for visualizing Micasense multispec point clouds. --- contrib/iterate_DTM.ipynb | 328 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 328 insertions(+) create mode 100644 contrib/iterate_DTM.ipynb diff --git a/contrib/iterate_DTM.ipynb b/contrib/iterate_DTM.ipynb new file mode 100644 index 000000000..133bdf6c4 --- /dev/null +++ b/contrib/iterate_DTM.ipynb @@ -0,0 +1,328 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + " Joshua W. von Nonn \\\n", + " Geographer \\\n", + " U.S. Geological Survey \\\n", + " Western Geographic Science Center \\\n", + " P.O. Box 158 \\\n", + " Moffett Field, CA 94035 \\\n", + " Email: jvonnonn@usgs.gov " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# SfM RGB Point cloud processing with PDAL\n", + "\n", + "The following code in this notebook will classify ground points and generate a DTM.\\\n", + "PDAL documentation: https://pdal.io/en/latest/\n", + "\n", + " Prior to running this notebook, the point cloud should have the outliers removed and thinned using PDAL's Poisson. I recommend a radius of 0.25 to start, but this will depend on your objective, microtopography, and the intended DTM resolution. " + ] + }, + { + "attachments": { + "image.png": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAANAAAACACAYAAAB6KueDAAAABHNCSVQICAgIfAhkiAAAABl0RVh0U29mdHdhcmUAZ25vbWUtc2NyZWVuc2hvdO8Dvz4AAAAtdEVYdENyZWF0aW9uIFRpbWUAV2VkIDA2IFNlcCAyMDIzIDA1OjE3OjAyIFBNIFBEVJ3XhhAAACAASURBVHic7Z15fFTlvf/fM8lk33eykJWQhZBACCFsQQRBQAQLFtyvhVttr7XVLvdebX+9t4u3tddrrdraal2qKAgCiuwxIQQIIQlJIPu+7+tkzyy/P0ImmUyWM5nBBHverxcvnTPnPOd5Ts73nOd55vl8PxJlb7UaERGRGSGd7QqIiNzJiAEkImIAYgCJiBiAGEAiIgYgBpCIiAGIASQiYgBiAImIGIAYQCIiBiAGkBH43btfcubyjdmuhsgsYDrbFfgmsCjIGy83h1k5d05xNacv5VDf0oEa8HZzZHv8UoJ9PQSX0T8wxJGEdLKLqhgYGsLTxYEd65ay0N9TcBnZRVUkphdQUddM/8AQLz2zG0c7a833L73zBZUNrTrHrV0awkP3rhB8nrmGGEBGYOuayFk9/4rFgXg4OyCVSknJKuL1g+f55Xd34GRvI+j4Iwnp5JXV8tSuu7C3teLclZu88elX/Pr7u7CzthBURl//IAt83Anx9eBYUqbO90/uWMugQjm6f98A/3fgLEsWzhfWyDmKXgGkVqv5/MJ1ruSU0N3Xj5OdNRtiF7F26UIAPjp5hZZOOc/uvUdzTHO7nJ+/eYQX92/H282J7r4BPj6VSkFFHQqlEg8ne3ZtXM6C+e6C6pCSVcSJ5CxeemY3EokEAIVCxY9f/ZjHt61mSYgvAPlldRxLzKCupQNHOyvWx4SxblmoppwhhZJnfvcP9u+M50pOKYVV9Zibynhi+yoWBflM21aAVz48TVFlAwA774pm08oIrbp29fRz4NQV8strkUgkLF7gw97NcViaywBISs8nKb2AmPAAkjMLUCrV3L08jHtXL9Zpd2d3H1Ip2FpZam1fvMBH6/N8DydSc0oor20WHEDltc1ELZxPgLcbABvjFpGUUUBDSwd21sLeZCsWBwForsd43J3ttT4npOXhZGdNiP88QeXPVfQKoIz8Si6kF/DU7vW4OdnR1NpJW1eP5vu4yCD+8MFJOrv7sLcZ/kOn3SzD280JbzcnAD5Puk5rh5znHt2MhUxGZX0LarXw9axRC305cCqV0pomgnyGgy63rBa1GhYFegNQ29zBm4cT2L0hhlB/L5raOnnv8xRsrCxZFuanVd6xpEx23hXNvp3xtHR2o1SpBLUV4LlHNgPDY6CJ+OBECvKePp579F6UCiXvn0jh4NmrPHHfas0+Te1dSCQSXnpmN8VVjbz60RmiFvoyz1X7hvvN28exs7bixf3bJ702A4MKktLzkUqleLk7CbiawwT5uFFQXk9XTz+2Vuak55bjYGvFfA9nwWXoS8r1QlZGLtA8BO9U9Aqg1g45draWLJjvjkQiwcHWSuv7AC9X3JzsuHazjA0rwgG4llvGmiXBmn1aOuT4eDjj5eoIgLODsKfkCDaW5oT4zSMzv1ITQBn55SwO9kEmMwHgVEo2seGBrF0aAoCroy3xy0K4nF2sE0Ax4f4sDR3e5u3mKLit09HW1cPNkhqef/RefG/diPfHL+XtYxfYc08sFrfeQjJTUzavXIRUKmWh3zwc7W2oamjRCaCpGBhU8KP//QiVSo21pTnP7tmAh7Pw43dtXM7Hp1P56aufIJVKcbSz4tmHNmnqaGyKKhtoaO0iLjLotpT/daJXAC0N9SMhLZdf/PkzQvw8CfJxY1lYACYmo0+RuIggUm+UsmFFOBX1LTS1y1keEaD5fs3SYN45doGaxjYCfdyICPTSa7AKEBPmz/ELmezeGINSqSanuJon71+r+b62uY365k5Ssoq0jpvopgrwcp1xW6eipV0OgJf7aFB6uzuhUqlp6ZRr3sj2NpZIpaOToRYyU3r6B3TK+/0P90x6LjOZCS/su5++/kEuZxfxj5OX+fHjW3CwERb0F68XUlzVwA/2bsTG2oLk9EJeP3iO//zOfVhbmAsqQx+SMwsJ9ffEacwkw52KXtPYro62/Op7u9i1IQZTUykHz1zlrSNfae0TuziQupZ2aps7SLtRRnigl1a/fclCX37z/d2sXbqQDnkvf/zkHAlXc/WqdFSIL929/ZTVNpNXXosECPP30tpn513R/OWFJ7T+/fKpnTplWZhPfIMIaeuU3OqWSpg64KQTdGH06NEOn0MiwcvVgSAfNx7bthpTUxOSMwsFHatUqvnsq3S2rFpMWIAX892deejeFfT0DZB2s0y/igigu6+frMJKVkctMHrZs4HevwOZm5kSGTyfb98Ty6P3reJGSQ2qW+MGAAcbK0L8PLmSXUR6XhkrF+u+pu1tLImLDGLfznjil4aQXVytVx0szWWEBniRmV9BZkElkQt9MTUdbYqnqyPFVY36Nk2H6do6FS5OdgDUNLZpttU0tiGRSHC2t9W7Lm1dPXR09wraV61WMzio0NqmVKpp7ejWebsNKRQMDSm1xiISiQQpEoaGlILK0IfLWSVYmpsRGew74zLmEnoFUGpOCVeyS2hs7aSprYv0vHK8XB21uiAAKyMXkJiej0KpIiJIe5boi+Tr3CypprWzm6qGVooq6zXdGX2ICfMnI7+CnMIqYsaNa+5dFUleeR1HEzOob+6kurGNc6m5nEsV/qYT2tbJcLKzJjzQi8Pn06hsaKWstpnjFzJZFu6vmYXTh//5+xe8/vF5ne0Hz1zlUlYxJdWNFFXU894XKTS1yYkM1p4ebm7v4oU3DvNlcrbWdgtzGUE+7py6lDM8Nmnp5NNz1xhQDBEW6CWoDIDuvgGqG9tobu8CoKFl+LoPjAlktVpNSlYhsRHCu8JzHb3GQObmMk5fyuGTs1eRSiDA2439D8Tr7BcVPB8zU1OiQ/213gww3KU5fD6d1s5uLM1lRAb7cv+6pXpXfPECH/5x8hIyExNC/LT/0N5ujjzz7Q18kXydr9LyMDM1xWeeM5tWLDJaW/PKannt43Oaz+V1zRxNzMDUVMrrP3sMgMe2rebAqSu88o9TSCQSwgO92LMpVu+2ToWdtSUJV3Np7epGIpEwz9mBp3evF/yzAMC+B+L5LCGdt48mMTCkwMPZnqd3rdeaVJmOzPxyDpxK1Xz+48dnAfjxY1sI8hmeHi+sqKepTc6qqIUTlnEnIrkdORG6evr599cO8dwjmzQzZSIi30SMuhJBpVIh7x3geFIGHk52YvCIfOMxagCV1bbwhw9O4u5sx5M7dLt2IiLfNG5LF05E5J8FUc4gImIAYgCJiBiAGEAiIgYgBpCIiAHcdkHd348nI5VKtZbwT8XfjiaRkVcBDGtMhB43lzmamEF1Qxs/2LtxtqsiYmTmnCJ1/8517N85HHjfFHw9nLG1FKbsNDZNbV18cuYq1Q1t9A8O4uJow/qYMNYsEb4aoKmtiy9TsimsqKenfwB3Rzs2rVxMTLi/4DKMIT2fi8y5APomMqI3mg1UKjXB893ZGBuOhYWM4spGDpxKxdrCXHC9KupaUChUPLhxOfa2VuQUV/POsQtYmpuyaNxax6kwVHo+F9ErgF4/eB5LCzPkPf1U1bVgb2vJ3s1xWk+RnOJqPj2XRoe8l8UL5qNUqbAwHx1qnb+ax5XsYpo75MhMTIgI9mH3xhij6U7UajX/+afDbFmzWOspeyolh/S8Mn7+rzuA4afqoXNpFFc1YGEmY8lCXx64exlmstFLcjQxg4r6VpaF+nLq0g26enoJ9PHgRw8NS9arGlo5dPYq1Y1tmEilzJ/nwnd2rNHIN05czOJEchYAYQFeE3bhUnNK+DIlm/auHlwcbdmxLpqoMXkCnn/lY+5aFkpeeR11ze3Mc3bgyR1rcXXUXtE9MKigt38AGysLZKYmmu0eLvZsdhmViPt7upKZX0FBRb3gAFq+KIDli0Y1XYHebtwoquZ6YZXgADKG9HwuovckQnpuOZvjFvG/z+9l7dKFvHU4kYHBIQA6e/r422dJRIf588K+7bg525JTXKV1fO/AAPfFL+HFfdv5wcP30NTWxYFTV4zTGoaX4i8N9SMjv1Jre2ZBBdFhwzfBkELJHz8+i7O9Df/x5H1878ENVNa3cvh8uk55dY3t3Cyp5flHNvP7H+4hfoy69u/Hk3F3tucX+3fw08e3EhXsg1I1+rv0tjVR/OWFJ3RyJYxQUd/C+ycusT4mlJ/vv5+YsAD++lkSTW1dWvul5Zby5P1refmHe7CwkHEsKUOnrKs3S/iPP31KXlntpNdGqVSTV1ZLXUsHvvNcJt1PCH0DQ9hYzaxbOjCoIPGa/tLzuYjeXbgFvh6E3FKQxkeHcDIlh8z8SuIig7iaU4q1pTn3xy9BIpGwfe0SHVHW9rVLtD5vWB7GhycvG9AEXZaF+5H4fh7dff3YWFrQ3C6nurGNfTuHlxel3ijB1ETKnk2xGh3MzvXR/OngOfZujtXSxgwqFDyxfY1GgjD2qd3a0c3W1VEaWbqHi3AZNcDFzEKCfT24KyYMGM7uk5lfzsXrhXzr7hjNfqsig3G5dY6YsABOXcrS84oM522oqGtGIpXwwF3LWGWAoO1Kdgl9g4PER4fodZyh0vO5iN4B5O5op/l/qVSKi4MNTR3DT8zmdjnzXB00N6BEIsHTVXtJfEF5HScuZlPb1EbfwPCby9jaEH9PVxztrMkqrGJ1VDCZ+RXMd3fWZIapaWynsbWLp3/7vs6xXb392FuPKmg9nO0m1e+sWxbC+19cJPVGCQFebiwN9WOeHkHU3C7X0UJ5ujvR1CbX2jZW+mxhZkpP76BOWWuXhmhyQEzEkzvW0tM/QFFFA1+mZOHj4cRCP/0z4hRXNXLw3FX274zXW5JtqPR8LqJ3ACnV2opMpVJ7KZ3JOMHZ2M8d8l7e+DSBLSsjefrB9VhbmJNdVMVbRxL1rca0RIf6kZFfweqoYDLyK4geJ7qbbEwyHospxmbfujuGuMgF5JfVkVVYyalL2Tz/2L34e06cZ2E8arUaIUlpxmeuUaP/8kVXR1tcscVvngvN7XJOXMzSO4DKapr486EEHt2ykvBbGZD0YUR6DsOZgH751lGSMwt1eiV3EnqPgWrHSJSHhpQ0dXTh6jD8VnJzsqO1Q/vp2dI+2p+vaWwFNWxeFaGZNKhr6pjwPOYymY6kWB+WhflTWNFARX0L1Q2tRIeNTrl6uTlQ3dDK4JBiihKE4eniwN3Lw3j+0Xtxd7Int3TyMch4XJ3sqGlq19pW19imM0EghP6BIVo7ugW1aSLJN4C8t5/Wju4J04yV1zXz+sHz7Nm8QutajsdQ6fmdht4BNCKPbmjt5NPzaUiRsPRWMsPYiEBaOrq5Xjg8gM8uqtK6QVwd7VAoVZTWNAFQ09ROUkb+hOfx83ShqLKeqoZWOnv6UCiE5SIYYb6HMy4ONnzweQp+ni6aMQRAXMQCzGSm/PVIEuV1LTS0dnL1Rin/+PKS4PKHhpQcPHOVkupGOuS93CytoaVTjo8eg+I1S4Ipqqgn8Vre8G8tF7Opa+lkdVTw9AePIy23lBfeOEx+eZ3W9qT0fM5euUlheR2lNU2cuXyDyzklRIXo5iT45EwqL7xxmL5bk0IjVDe28dqBs6yMDGKeqwPVjW235NtynTIMlZ7faejdhVu+KJD88jqOX8jA1cGWp3av1+QPs7O24F+/tY5DZ9M4ePoqfl6uWjkR3J3teXhLHO8ev4hKpcbZwYYNsYs4mqg7+xUbEUhxdSP/+4/TDAwO8fTu9Xpf7OgwP06l5LBrzIAcQCYz4dmH7uFIwjX+eOAMarUadyd7VuqRp0wildDZ08c7R5OR9/XhaGvFfWuWaOo4kvl0LE/95j0AfvrEVgK8XPH3dOXRras4dTmbwwnXcHGwZd+OtXpPRkyFg60151JvcvpyDkMKJc4ONuxYt5SNt/L2CaGgvJ6+gSHOX83j/NU8zfaF/p6aKf3pMIb0fC6ilx7o9YPncXO048F7lt/OOomI3DGIi0lFRAxADCAREQMQJd0iIgYgvoFERAxADCAREQMQA0hExADEABIRMQAxgOYAQwolT/3mPcpqm2e7KiJ6IipS5wBSiZQVi4OwnaG+xlDqmzs5fiGTitpmOrp7+Zfta4iNCBR8/OCQgi+Sr5NTXE1bZw/2NpbERgSyZVXUhCvtc4qrefNQAtFhfuzfuW5GdTZGGcZADKA5gImJZFaTp/QPDeHqaMPycH/++lmS3sfLe/spq2nmnhURuDvb0djayaFzaQwOKbR0TcP79nH4/DW9nB90z2d4GcZCrwB660giZjJTFEoleaW1WFua89C9cYQFDNuL1DZ38Ku/HuN3P/y2RlPzyelU2uU9PL37bk5czKKgvA55bz9SqYS7loVx8mI2TvbW/GDvPViYyziXmktqTjHBvvO4lluOVCphy+rFGoft9NxyPjx1mZd/uEcjXVar1bz4xhHujg1j/S1x2lS8czwZ1Gq+MyZ/d0l1I698eJr/+cG3NdbuCWl5JF7Lo7OnD09nB7519zKCx0gApnP6FuJIPrI+DkbXyI3FWE7fKpWKzu4+LMzNdPRN/p4u+HvOXKHqbG/DTx7fovkc5ONOdUMbWYVVOgH0wYnLbI6LID2/YsbnM0YZxkLvMVBGfjlrooJ55fmHWLE4iA++uCTYtQ2GNUHff3ADJhIpF9ILeO7RzXT39XNjjEtdbXMHNlYWvPyjb/P0g+s5kpCuWcEduXA+UiTkFI1KxUuqm2iX9xITHqBzvomICfMnp7iGIcWoXOJ6QSULfOdpgic5s4Dzqbns2RTLL/bvYGXUAl4/lEBrR7dOeceSMlkZGcTLz+7hh49swvaWQGysI/kv9u9g08oIHanAX154gj/97NFJ6/rBiRTau7p57tF7+cGejVTWt3Dw7FWtfcY6fe9/IJ7PL2RS39ypvU+bnP/406d8ceG6oGtkKL0Dgzpd0ovXC1Eolaw0QA1rjDKMid4BFDx/HiH+nkgkElZEBNHR3UtHd5/g4308nHFzssPX04UAbzfcnOzwcnOipXP0xjSTmXJP3CIkEgn+nq6EB3px8fqwYbDM1IToMH9Sb4xKxdNulrEo0EvwGCIswAupBE3+ALVaTWZBpZbT3YnkbB5YH82iIJ9hl+/oEPw8XUjPK9cpb8Tp28Jchrebo8aVe6wjubODDUtD/fRK4zTi9L1rw3J8PZwJ8Hbj/vilpN0spX9gVHIwmdP3bFHf0sn1wko2rxp9Cza1dfHFhSwe2bJyxuUaowxjo/cYSEtebD58eG//gGB570i3S2ZiorFLlJmaMKQYFVY52llpZZZxc7KjrGZ0hiouMog/fHAKeW8/luZmZOSX89g24WMIUxMpkQt9Sc+vIDJ4PuV1LXR29xG1cFgjI+/tp6unj3eOJw9398bg7iTc6dtQR3JjOn17uNjzlxeeEHzumdLZ08efPz3PhthwTSYetVrNu8cvsm1tFM4zzMBjjDJuB3oH0Hh5MYy6Sk8kT1YJtJweu9f4LuF42XiAlytuTrak55XjZGeNRCJhUZB+EuOYMD/+dvQCQwolmfnlhPjNw8ZSW779s3/ZKkiePZnT94gjeV5pLblltfzxk3N8a/0y7o4VqMX5Gp2+jYG8t49XPzxNeIA398eP2nYODikor2umuqmVQ+eGu58jf9N/K/qA137yyLTes8Yo43Zg1Fk4K3MzAAYGhuDWJEJbVw8mUv2ShrR19dDbP4iVxXB5tc3tuIyTOcdFBJF2oxQHO2tiwgMwNdHv4oX4eWEikZJXVsv1giq2rYnSfGdrZYGdtQXFlY2C8xtMxogjeVxkEAfPWJBdXC04gMY6fY9MPMzU6VupVNMh78HCQjbjHHxTldHd18+rH54hyMdDRy9mJjPll9/dqbXtw5OXsDAzY9eGGJ0bv62rB6lUopVsRN8yvi6MelZ7G0scbK24ljs8Tiirbaawsl7vclQq+PhMKg0tnaRkFVFYUa8jc45dHEhlQys5xVWsWCz8N4sRTEwkRIX4cjwxk87uXiJDtNWuW9dEcfJSNpeyimlul1NS3cihs2nkl9VNUqIuhjqSG9PpeyqHbYVCpZFpw/ANXN3YRldPv6Ay+gaGePXDM1hbWbAmOpiapnaqG9uoaRouTyKR4OFir/XPTCbD3Nx0QvXtRLJwfcv4ujDqG0gikfDYtlUcOHWFpPR8gn3nER3qT9+4/vh0eLk6YmdlwUvvnsBcZsqD9yzXOD2P4GBjRYifJ62dcvxmmCQwJsyPS1lFRAT5aN6eI8RHh6BSqThz5QYHTl/G1tKSwPluOm/CqZjOkfxoYgZnLt/QfP79e18C2hmDvg6n75YOOb95+3PN5+NJmRxPymR7/FK2jJsOn4imti5N7ovfvv2FZrtUKuHN/3jcqHWda8w5PdC51Fyu3ijlxf3bp933V389zrIwf53fPEREvi7uyJUI3X39XC+opLGtkzg9EoGIiBibOzKAfvrqQazMzXn43pU42N65WS1F7nzmXBdOROROQpQziIgYgBhAIiIGIAaQiIgBiAEkImIAd0QA/ftrh7iSXTLb1ZhVjiZm8NrH52a7GiLjuCOnsf8ZmW2nb0NdugFullRzJCGd5g45Hk4O7Nm8QmeFyVTMRadvMYDuEGbT6dsYLt2NrZ385XAi8dGhxEUGkXgtj9cPnuNX33tAY8oshLnm9G1Ul26VSsVHp65QWFFPh7wXW2tLVkUuYOuaSI0MQojr9JcXs0m8lo9KpWLjikU69RDi9H0tt5yTKVk0d8ixtjAnIsiHR7YKE2I1tXXxiz9/xn89tVNjCwnw6oGzeDjba9aiCXH6/t27XxIW6Elv/yBXb5SiVKm4OzZc48o2XT0Ndfoekdk/vCWO81dz6ezuIzzQiye2rUEmM9Eqp6d/gMFBBQ62VlqyFWO4dF/OKcHRzopdG5YhkUh4aPNKsgqquHqjnA2x08vwYW46fev9BkrPLefZvRtZ6DePpPR83jqcyG+f2YW5mQyVejjDzGNbV+HsYEtTWxfvHEvGzsZCy78zLbeUH+zdhL2NJW9+msCxpAxNZpWswipOXc7m8a2r8fJw4rOEdLp6tBWvI07fnq4O9A0OcfDMVQ6cuqIpo0Pey7ufX2Tv5lgignzo6umjoFz4qnA3JzvmuzuTmV+pWWfX3TdAUWWDRvYw4vS9KNCbXU/ex8Cggk9Op3L4fDoP3btCq7zkjCLWLQvhv7/3AAqFirrmdsH13LYmim1rojiamEF1QxvjGXH6fnBjDGEBXqTnVfDXz5L45Xd34OY06md7o7iG/3zyPvoGhnjp719w5UaxjqfqgVNXyMir4JUfP6SzuHY8+rp0V9S14O/lpglMExMJ8z2dqaibWSqvgUEFSemz7/St9yTCiEu3RCIhPjoEqVRC5i1LeVMTKQ9viSPYbx7ODjaEBniyLMyPvHESgBHXaZmpCTFhAVTVj8qPU64XsSTEj5hFAXi6OLB30wpUKu3FEtvXLiFq4fxhabiHMxuWh2nJDNrlPahUKhYvmI+DrRXzPZy5J073TTYV0WHDHqsjZBdWYmdtQaD3sD5orNO3h7M9vvOc2bk+mis3inXyHrg727J1TSQ2lhY42FppkrAYo55jnb7dne3ZuiaSec52XLxeqLXf5lWLMTeT4WBrRWiAFxX1rXqdZywzcenu7unH1sqC/LI6nvvDASrrW7G1skDe2z/9wWMYGFTwvZfe59mXP+Rs6s1Zd/o2qks3DGeJuZBRQEtHtyZpR2iAtox5Ktfp5o4u4nxGF4g6O9hoHPBGmM7pe777cL6F/37rGKEBngR4uxITHqBX3rVlYf4cTcygqa0LNyc7MvIriQ711zxB9XH69p9E8m2Megp2+rYfc81lprTLdX1M9+9cx/6dOpu1mLlLtxqlUoWlpRlO9tbIZCY6SmMhzDWnb6O6dGcVVvHZVxns2xlPiN88zGSmHDqbRm2LtpHwVK7TEiSYmGj3zU0k+jl9m5hI+PGjmympbqKosoGk9HzOXcnl/313h04wToazgw3+ni5k5FcQHx1CYWUd98VHae0j1OnbchLJtzHqOVOnb2bg9G2IS7eNtSXyvn785rnw4v77AejuH9A7meRcc/o2qkt3aU0jQT7uLF7goxlI17VM7MI9GS6OtjSPcfbu6R+gd4wgT6jTt1QqJdjXg21ro/jpE9tol/dQ1aBftyU6LIDM/AqyCqtwsLHWkncby+nb0Hoa0+n7drp0+85z1koMo1SqqaxrnlAMeSc5fRvVpdvD2Z66pnZNvzY9t5yiyga9yl+7dCEZ+RU0tXWhVqs5eTFH61kpxOm7vK6ZM5dvUNPURntXD5ezipGZmmjNqAkhOtSXmsY2zqfeJHpMyiswjtO3MeppTKfv2+nSvSpyAfKePg6fT6O2uZ0Dp68AELtYN5ffneT0bVSX7rjFQVQ1tPHrvx1HKpUS4O3K6qhgGse8UaZj8QIfNiwP5+X3T2Jhbkp0WIBW/1aI07eFzIz8inrOpt5kUKHA85YjtL2N8N8bABztrAn0caekupHH71+j9Z0xnL6nq+dccfo2hku3u7M93911F0cS0klKL8DD2Z7vP7hBr9+A5qLTt+jSLSJiAHfEWjgRkbmKGEAiIgYgSrpFRAxAfAOJiBiAGEAiIgYgBpCIiAGIASQiYgBzNoBE52r9EK/X7DBnFanGcK7+9d8+JzYikI0rBPrx3MHMttP3wOAQH5+5SlZhJVKJhOXhAezeGDuhS/d0qFQqXn7/JOV1Lfzyuzs1qyqGFEoOnL5CRW0zDS2drIxawKNbV+lVtqGO5OOZs2+gEefqmSyK/Gdktq/XJ2evUlBWx/d2382+nfFcyy3n8+TMGZV18lIOpqYmOttVKjWmUhM2xi7Ce4YiuhFHcmOtpjGqpPujk1do6ZTz7N7R9VHN7XJ+/uYRXty/nYq6Fk4kZ/HSM7s1y+sVChU/fvVjHt+2miW3FqVO51w9nfv1z988olnoWJPQxpGEawA8vXs9kcHzBUnPp0OI0/fAoIIjCdfIzK9EqVIS6OPO3k0rcHYYlR/nldXyQ7fCwAAACbVJREFU+sHzPPfIZg6euUp9Swf2tlb85JbGxRhO3+V1LRw6c5XqplasLcxZuzSELasXa9o6nfv6CJM5ffcPDJF2s5THtq7W3Av3ro7kZEoW98cv0cv8qryumYz8cvbes4JXPjqj9Z25mSkPb4kDIPVmqeAyx2KoI/l4jCrpHvYuPUlnd59mQWTazTK83ZzwdnPCwdaaA6dSKa1pIshn+AbILatFrYZFY/Qlf3nhiQkXUo4w1v3aQiajsr5Fawn+r773LWDyLpxQ6flUxIT5886xZIYUSo2f63in7w++TKGru4/v79mAlbmMs6k3efNQAi/su0/rplKr4VhiBg9tWYmniwOlNY2Y3Pp+urZOd70GhxS8eTCBiAVePH7fauqa23n/RAqOdtZaC18z8sv5twc3sG9HPF+mZPPBF5f47TO7tOrZ1Cbnl28dZX1MmNYTvLqxDaVSrSUcDPByobd/kKZ2uWDF6MCggvc/T+GRrauQ6ulqOFsYVdI97F1qx7Wbow7a13LLiLvlIGdjaU6I3zzN/jD8h1sc7KOT4GIqDHW/Fio9n4rpnL4bWjvJzKvgOzvi8fd0wd3Znr2b4mhs66RynN5HrVZz35oo/D1dMDczJSxg1HHc0LZeL6hkUKFgz6YVeLjYszTUj7iIIC6Mk4AY4r7e3TcsX7G1tuSVD0/z18+SNPkS5D3CJduHz18jNMBTr1RXs43RJd1xEUGk3ihlw4pwKupbaGqXszxiVPMRE+bP8QuZ7N4Yg1KpJqe4mifvX6tXHQx1vwZh0vOpmM7pu66pHTXDSSHH09rereO96uc18U1jaFubO+S4OthqZQry9nAibcxDDoS5r0/q9H3rhahUKYcl+DIZSpV+K8RultZQUFHLi/t26HXcbGNUSTcMe5cev5BJbXMHaTfKCA/00tJ8RIX48tGpy5TVNtPTN4AECPPX7mtPh6Hu10Kl59MxndO3TGbCaz95ZNpxlVQqwdxs4j+FoW2dueRbuNP3yNumu2eAx7etBqCkeljwaGstbFawoLye1s4env+/A7dOPvyfX799nM0rF7NtbdTkB88iegfQRJLudQ6hmm0j3qVXsotIzytj7+Y4reMtzWWEBniRmV9BT/8gkQt9MTXVfzJQiPu1zFSKUqnUOXas9HwEfaXnMLXTt6erI0NDSiobWmfs4TqCIU7fbo52NLXLGRxSaN5CNQ0zk3xP5tLt4+6EiYmE0pomzZRzaXUjVuZmuI07T//AED19A9haW2i9FTetjNBS0dY0tfP20SSe3r0eXz0H/cZwJBeKUSXdI6yMXEBiej4KpYqICRLvxYT5k5FfQU5hlWbMoA9C3a9dnOzIK6+jq6efIYVSM/g2hvQcpnb6HhlvvHs8mbyyWprb5eSW1vCXw1/RPzA0Rakza+tkRC30xczUhE9uuZ5n5ldw5UYJa6IX6tVWmNyl28JcxvJFgZxIzqKosoH8sjrOpt5kTfRCnRm4tNxSXnjjMPnl2uNNWysLLQfukSxCzva22IxJaVzX0kF1YxsDtwJxImm5MRzJhWJUSfcIUcHzMTM1JTrUf8K3y+IFPvzj5CVkJiaE+Gl334Q4V0/nfj3C1tWRfHDiEi+88SlDQ0rNNLYxpOcjTOX0/fi2VRxNzOS9z1PoHRjAyc6GsACvCX/jmAxDnb7NzUx5evd6Pj2Xxq/fPo6VhRl3Lw9jVeQCvds6FXvuieXjM1d589MEJAzfJ/fdhkw5r354WnOzVza0klVYpZe03FBH8vHcFkl3V08///7aIZ57ZJNmulpE5JuIUZfyqFQq5L0DHE/KwMPJTgwekW88Rg2gstoW/vDBSdyd7XhyzC/0IiLfVERJt4iIAczZxaQiIncCYgCJiBiAGEAiIgYgBpCIiAH80wfQ7979UuuHyG8y/0xt/bqYs5Lur0uOvSjIGy83h9t6jrnCbLbVGHJsQ126BwaHePbljyb8bt/OdSybwbKyORtAXxdb10TOdhW+NmazrWPl2Inp+dMfMAmGuHSbyUx5Yd92rW2F5fUcT84kPFA/RcAIegfQ0cQMKupbWRbqy6lLN+jq6SXQx4MfPXSPIKm0oXJswChS6Vc+PK1ZQLrzrmg2rYzQaqdKpeJoYgapOWX0DQ4Q4OXGns0r8HQR9gRPySoSJF/PL6vjWGIGdS0dONpZsT4mjHXLQrXKMtTp29C26uP03dndh1SKjm2JMeTYhrp0SyQSfMblUjiRfJ2lIX5aEnV9mNEbqK6xnZsyU55/ZDOWlmYU3FJyCpFKGyrHBuNIpZ97ZDMwfHNOREJaPpezS/iX7atxtrfjWFIGfz6UwH89tVOQxj9qoe+08vXa5g7ePJzA7g0xhPp70dTWyXufp2BjZanTnTDE6dtYbRXi9P2bt49jZ23Fi/u3jz+NUTGGS3eHvJcbJTX86OHNM67HjCYRBhUKnti+BmcHG6zMzVga6gcIk0obKlE2llR6OpIzC1gfE8aiIB/mudrz8JaVtHZ2ayTc0yFEvn4qJZvY8EDWLg3B1dGW8EBv4peFcDm7WKe82+n0LbStxnT6ninGdOlOuV6Ei4OtQQZdM3oDeTjbTfrKm04qbahE2VhS6alQqVS0dnbj7e6o2WZnbYGdlSVNE9gaTsZ08vXa5jbqmztJySrSOm6iG+J2OX3r01YhTt+//+EeQeedKcZy6VapVFzKKmbdMmFJZCZjRgFkMYnKT4hU2lCJMhhHKi0MwzLDCJGvTzQmmYjb6fQ9zPRtNYbTt6EYy6X7ZmktnT19rNDDknMijPo7kFCX7hGJ8r6d8cQvDSG7uFpnn8nk2GOl0rcLqVSKs701NWPk6109/XT19uklhR4rX88sqNSRr3u6OlJc1WiU+s7U6dtYbR1BH4dtY5QxkUu3UqmmtaObnjHu7uO5kFFARJAX9tb6+eaOx6gBJEQqbbAc20hS6elYsySEr67lcbOkmvrmTj46eRlHW2vC/L2nP3gMU8nX710VSV55HUcTM6hv7tTI5c+l5gou3zhO38ZpK0zusA3C5NhTlSHUpXsqSTdwa3xXx6oZOJmPx6i/AwmRShsqxwbDpdJ5ZbW89vE5zefyumaOJmZgairl9Z89BsDdy8Pp6unlvc8v0T84iJ+nC08/uF7vXM9Tyde93Rx55tsb+CL5Ol+l5WFmaorPPGc2rRA+CTCd0/fX2dbpMFSObSyX7pSsIuytLVk0w99+xiLqgUREDOCffi2ciIghiAEkImIAYgCJiBiAGEAiIgYgBpCIiAGIASQiYgBiAImIGIAYQCIiBiAGkIiIAYgBJCJiAGIAiYgYgBhAIiIG8P8B4JFr71a2ZA8AAAAASUVORK5CYII=" + } + }, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "![image.png](attachment:image.png)" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "sys version:3.8.17\n", + "pdal version:3.2.3\n", + "numpy version:1.24.4\n", + "pandas version:2.0.3\n", + "pyvista version:0.41.1\n", + "rasterio version:1.3.7\n" + ] + } + ], + "source": [ + "import sys\n", + "import pdal\n", + "import numpy as np\n", + "import pandas as pd\n", + "import pyvista as pv\n", + "import rasterio as rio\n", + "\n", + "print(f\"sys version:{sys.version[:6]}\")\n", + "print(f\"pdal version:{pdal.__version__}\")\n", + "print(f\"numpy version:{np.__version__}\")\n", + "print(f\"pandas version:{pd.__version__}\")\n", + "print(f\"pyvista version:{pv.__version__}\")\n", + "print(f\"rasterio version:{rio.__version__}\")" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "#normalizing band color values for pyvist\n", + "def normalizing(dflist):\n", + " dflist_norm = (dflist - dflist.min()) / (dflist.max() - dflist.min())\n", + " return(dflist_norm)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "#enter absolute path of laz or las file, example: fp = '/home/user/Downloads/points.laz'\n", + "fp = '/media/grendel/AE10528B10525B03/UAS_2021_Sonoma_Fires/GLA3_2021/WebODM/1cm_products/Glass_ultra_croped_and_thinned_olr_rmvd_byhand.las'" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "'Pipeline selected 633600 points'" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#reading las file into pdal pipeline\n", + "p = pdal.Reader.las(fp).pipeline()\n", + "n_pts = p.execute()\n", + "\n", + "f'Pipeline selected {n_pts} points'" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "#uncomment to thin point cloud to 25cm radius\n", + "#p = pdal.Filter.sample(radius=0.25).pipeline(p.arrays[0])\n", + "#p.execute()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### To iterate, set the parameters for the algorithm of your choice then skip to \"Setting viz parameters\"" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "633600" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#smrf parameters ---- settings for complex topography\n", + "\n", + "cell = 1\n", + "window = 13\n", + "scalar = 1.5\n", + "slope = 0.35\n", + "threshold = 0.5\n", + "\n", + "ground_test = pdal.Filter.smrf(cell=cell,window=window,scalar=scalar,slope=slope,threshold=threshold).pipeline(p.arrays[0])\n", + "ground_test.execute()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#PMF ----- aggressive parameters for flat topography, i.e. meadow\n", + "\n", + "cell_size = 1\n", + "max_w = 33\n", + "slope = .05 #\n", + "initial_dist = 0.05\n", + "max_dist = 2\n", + "\n", + "ground_test = pdal.Filter.pmf(cell_size=cell_size,max_window_size=max_w,slope=slope,\n", + " initial_distance=initial_dist,max_distance=max_dist).pipeline(p.arrays[0])\n", + "ground_test.execute()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#cloth ---- default settings\n", + "\n", + "resolution = 1\n", + "threshold = 0.5\n", + "step = 0.65\n", + "rigidness = 3\n", + "\n", + "ground_test = pdal.Filter.csf(resolution=resolution, threshold=threshold, step=step,\n", + " rigidness=rigidness).pipeline(p.arrays[0])\n", + "ground_test.execute()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Setting viz parameters" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "'535024 ground points found.'" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "arr = ground_test.arrays[0]\n", + "cols = [col for col, _ in arr.dtype.descr]\n", + "\n", + "df = pd.DataFrame({col: arr[col] for col in cols})\n", + "\n", + "bands = (\"Red\",\"Green\",\"Blue\")\n", + "for band in bands:\n", + " df[band] = normalizing(df[band])\n", + "\n", + "df.loc[df['Classification'] == 2, bands] = (0.98,0.78,0.44) #colorizing ground points\n", + "\n", + "pc_norm = np.dstack((df.X, df.Y, df.Z,))\n", + "pc_colors = np.dstack((df.Red, df.Green, df.Blue))\n", + "\n", + "pc_norm = pc_norm.reshape(pc_norm.shape[1], pc_norm.shape[2])\n", + "pc_colors = pc_colors.reshape(pc_colors.shape[1], pc_colors.shape[2])\n", + "\n", + "pcn = pv.PolyData(pc_norm)\n", + "pcn['colors'] = pc_colors\n", + "\n", + "f'{len(df[df.Classification == 2])} ground points found.'" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plot will pop-out window" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "pcn.plot(scalars='colors', point_size=2, rgb=True, notebook=False, full_screen=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Write out DTM\n", + "https://pdal.io/en/2.6.0/stages/writers.gdal.html\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + " The DTM will likely have no-data holes. I recommend the GDAL fill no-data tool in QGIS to resolve these. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#selecting only ground points to generate DTM\n", + "#resolution here is set to 50cm\n", + "#DTM will be written to same location as this notebook\n", + "\n", + "DTM = (pdal.Filter.range(limits=\"Classification[2:2]\").pipeline(ground_test.arrays[0]) |\n", + " pdal.Writer.gdal(resolution = 0.5, radius = 1, output_type='min', filename=\"DTM_test.tif\").pipeline())\n", + "DTM.execute()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "pdal38", + "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.17" + }, + "orig_nbformat": 4 + }, + "nbformat": 4, + "nbformat_minor": 2 +}