Merge pull request #7 from SnowEx/canopy-file-filestructure

Looks really clean Zach!
cryogars · Sep 14, 2022 · d635d7c · d635d7c
2 parents 6badec3 + 0657084
commit d635d7c
Show file tree

Hide file tree

Showing 16 changed files with 479 additions and 166 deletions.
diff --git a/LICENSE b/LICENSE
@@ -1,6 +1,6 @@
 MIT License
 
-Copyright (c) 2022 Brent Wilder
+Copyright (c) 2022 Brent Wilder, Zach Keskinen
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal

diff --git a/README.md b/README.md
@@ -16,17 +16,30 @@ $ conda activate iceroad
 #### Running the code
 from the ice-road-copters directory for example one can run:
 ```
-$ python ./scripts/ice-road-pipeline.py /Users/brent/Documents/MCS/mcs0407 -e /Users/brent/Documents/MCS/QSI_0.5M_PCDEM_USIDMC_20210917_20210917_WGS84_proj.tif -a /Users/brent/Code/ice-road-copters/ASP/bin -s /Users/brent/Code/ice-road-copters/transform_area/hwy_21/hwy_21_utm_edit_v2.shp
+$ python scripts/ice-road-pipeline.py <path-to-directory-of-laz-files> -e <path-to-user-supplied-reference-dem> -a <path-to-ASP-directory> -s <path-to-road-shapefile-to-clip-to>
+```
+NOTE: this code assumes you are using a reference DEM (and heli lidar data) in ellipsoid. If your reference DEM is in Geoid you must set the `-g  True` flag. Also, this code is assuming you would like a buffer of 2.5 meters on either side of the centerline of the road (see buffer_meter variable in laz_align.py). TODO: turn this buffer_meter variable into a flag as well.
+
+#### Flags
+
+```
+-e user_dem      Path to user specifed DEM [one will be downloaded from py3dep if you don't supply one]
+-d debug         turns on debugging logging [Options: True or False]
+-a asp_dir       Directory with ASP binary files [Can be either ASP or ASP/bin directory]
+-s shp_fp        Shapefile to align with [road shapefile to use to tie reference DEM to your point cloud]
+-g geoid         Is the reference DEM in geoid? Will be auto set to True if you don't supply a DEM [Default: False]
 ```
-NOTE: that this code assumes you are using a reference DEM (and heli lidar data) in ellipsoid. If your reference DEM is in Geoid or you are using py3dep functionality, you must change `dem_is_geoid=False` to `dem_is_geoid=True`. **Additionally, HWY-21 is hard coded for now into the align function but can be easily changed in the future if this code is used for say Bogus Basin Road (see `transform_area` and `buffer_meters` parameters).**
 
 
 ###  Additional information :books:
-The goal of this program is to utilize existing USGS 3DEP 1m topography data (via [py3dep](https://github.com/hyriver/py3dep)) and Ames Stereo Pipeline ([ASP](https://github.com/NeoGeographyToolkit/StereoPipeline)) software to accurately align snow-on airborne lidar point clouds to real world coordinates without the use of ground control points.
-![heli_bsu](./docs/heli.png). We also provide an option for a user specified DEM (snow-off).
+The goal of this program is to utilize existing USGS 3DEP 1m topography data (via [py3dep](https://github.com/hyriver/py3dep)) and Ames Stereo Pipeline ([ASP](https://github.com/NeoGeographyToolkit/StereoPipeline)) software to accurately align snow-on airborne lidar point clouds to real world coordinates without the use of ground control points. We also provide an option for a user specified DEM (snow-off).
+
+![heli_bsu](./docs/heli.png) 
+
+For example, we used prior knowledge that HWY-21 running through our study site is kept snow-free for a majority of the year, thus making excellent virtual ground control points for post-processing in ASP (using `pc_align`).
 
-For our case, we used prior knowledge that HWY-21 running through our study site is kept snow-free for a majority of the year, thus making excellent virtual ground control points for post-processing in ASP.
 ![roads](./docs/roads.png)
 
-After running pdal processing and ASP post-processing, we are able to generate accurate snow depth maps with sub cm road differences
-![snow](./docs/snow.jpeg)
+After running pdal processing and ASP post-processing, we are able to generate accurate snow depth maps with sub cm road differences.
+
+![snow](./docs/snow.jpeg)
diff --git a/notebooks/align_check.ipynb b/notebooks/align_check.ipynb
@@ -51,7 +51,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3.8.13 ('iceroad')",
+   "display_name": "Python 3.8.13 ('heli')",
    "language": "python",
    "name": "python3"
   },
@@ -69,7 +69,7 @@
   },
   "vscode": {
    "interpreter": {
-    "hash": "14aac5828947b75d97bfce66d0ff9606be49ebe3317095af148ca98ecf57d790"
+    "hash": "b38bf96d929ed892a854a240f79fc97050674deef20950d0719a1a6703ca0f3a"
    }
   }
  },

diff --git a/notebooks/buf_laz_check.ipynb b/notebooks/buf_laz_check.ipynb
@@ -40,7 +40,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "gdf.to_file('testing/test.shp')"
+    "buf = gpd.read_file('/Users/zachkeskinen/Documents/ice-road-copters/test/mcs-data/ice-road/results/buffered_area.shp')"
    ]
   },
   {
@@ -49,7 +49,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "gdf.geometry"
+    "with laspy.open('/Users/zachkeskinen/Documents/ice-road-copters/test/mcs-data/ice-road/results/clipped_PC.laz') as f:\n",
+    "    hdr = f.header"
    ]
   },
   {
@@ -58,7 +59,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "buf = gpd.read_file('/Users/zachkeskinen/Documents/ice-road-copters/data/results/buffered_area.shp')"
+    "buf['geometry'].plot(linewidth = 100, color = 'black')"
    ]
   },
   {
@@ -67,8 +68,51 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "with laspy.open('/Users/zachkeskinen/Documents/ice-road-copters/data/20220317-220317_185608_VQ-580.laz') as f:\n",
-    "    hdr = f.header"
+    "import rioxarray as rxa\n",
+    "dem = '/Users/zachkeskinen/Documents/ice-road-copters/test/mcs-data/ice-road/results/ref_PC.tif'\n",
+    "d = rxa.open_rasterio(dem)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "d.plot()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "buf.plot()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import fiona\n",
+    "from shapely.geometry import mapping\n",
+    "import os\n",
+    "# Define a polygon feature geometry with one attribute\n",
+    "schema = {\n",
+    "    'geometry': 'Polygon',\n",
+    "    'properties': {'id': 'int'},\n",
+    "}\n",
+    "# Write a new Shapefile\n",
+    "with fiona.open(os.path.expanduser('~/Downloads/laztest.shp'), 'w', 'ESRI Shapefile', schema) as c:\n",
+    "    ## If there are multiple geometries, put the \"for\" loop here\n",
+    "    c.write({\n",
+    "        'geometry': mapping(laz_box),\n",
+    "        'properties': {'id': 123},\n",
+    "    })\n",
+    "# type(laz_box)"
    ]
   },
   {
@@ -77,20 +121,30 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "with laspy.open('/Users/zachkeskinen/Documents/ice-road-copters/data/results/clipped_PC.laz') as f:\n",
+    "with laspy.open('/Users/zachkeskinen/Documents/ice-road-copters/test/mcs-data/20220317-220317_172628_VQ-580.laz') as f:\n",
     "    hdr = f.header\n",
     "f, ax = plt.subplots()\n",
     "laz_box = box(hdr.x_min, hdr.y_min, hdr.x_max, hdr.y_max)\n",
     "x,y = laz_box.exterior.xy\n",
     "ax.plot(x,y)\n",
     "buf.plot(ax = ax, color = 'red')\n",
-    "ctx.add_basemap(ax = ax, crs = buf.crs)"
+    "# d.plot(ax = ax)\n",
+    "# ctx.add_basemap(ax = ax, crs = buf.crs)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "laz_box"
    ]
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3.8.13 ('heli')",
+   "display_name": "Python 3.6.8 64-bit",
    "language": "python",
    "name": "python3"
   },
@@ -104,11 +158,11 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.13"
+   "version": "3.6.8"
   },
   "vscode": {
    "interpreter": {
-    "hash": "b38bf96d929ed892a854a240f79fc97050674deef20950d0719a1a6703ca0f3a"
+    "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
    }
   }
  },

diff --git a/notebooks/pipeline.ipynb b/notebooks/pipeline.ipynb
@@ -326,15 +326,27 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "PC = laspy.read('/Users/zachkeskinen/Documents/ice-road-copters/test/09_EXPORT/crop.laz')\n",
-    "\n",
-    "points = np.vstack((PC.X, PC.Y, PC.Z)).transpose()\n",
-    "\n",
-    "cloud = o3d.geometry.PointCloud()\n",
-    "cloud.points = o3d.utility.Vector3dVector(points)\n",
-    "#pcd.colors = o3d.utility.Vector3dVector(colors/65535)\n",
-    "#pcd.normals = o3d.utility.Vector3dVector(normals)0\n",
-    "o3d.visualization.draw_geometries([cloud])"
+    "import numpy as np\n",
+    "import pdal\n",
+    "import json\n",
+    "import laspy\n",
+    "import open3d as o3d\n",
+    "import os\n",
+    "import matplotlib.pyplot as plt\n",
+    "import rioxarray as rxa\n",
+    "from os.path import join, basename, exists, dirname, abspath\n",
+    "from glob import glob\n",
+    "import geopandas as gpd\n",
+    "import fiona\n",
+    "from zipfile import ZipFile\n",
+    "import shlex\n",
+    "import subprocess\n",
+    "import py3dep\n",
+    "from shapely.geometry import box\n",
+    "import pyproj\n",
+    "from shapely.geometry import Point\n",
+    "from shapely.ops import transform\n",
+    "import contextily as cx"
    ]
   },
   {
@@ -343,7 +355,17 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "abspath('~/test/09_EXPORT/crop.laz')"
+    "import laspy\n",
+    "import numpy as np\n",
+    "PC = laspy.read('/Users/zachkeskinen/Documents/ice-road-copters/test/mcs-subset/ice-road/results/clipped_PC.laz')\n",
+    "\n",
+    "points = np.vstack((PC.X, PC.Y, PC.Z)).transpose()\n",
+    "\n",
+    "cloud = o3d.geometry.PointCloud()\n",
+    "cloud.points = o3d.utility.Vector3dVector(points)\n",
+    "#pcd.colors = o3d.utility.Vector3dVector(colors/65535)\n",
+    "#pcd.normals = o3d.utility.Vector3dVector(normals)0\n",
+    "o3d.visualization.draw_geometries([cloud])"
    ]
   },
   {

diff --git a/scripts/full-ice-road-run-SNOWDATA.sh b/scripts/full-ice-road-run-SNOWDATA.sh
@@ -0,0 +1,48 @@
+#!/bin/bash
+
+ASP_DIR='/SNOWDATA/IDALS/ASP/bin/'
+MCS_SHP='/SNOWDATA/IDALS/ice-road-copters/transform_area/hwy_21/hwy_21_utm_edit_v2.shp'
+DC_SHP='/SNOWDATA/IDALS/ice-road-copters/transform_areaBogusBasin/Roads_DryCreek.shp'
+MCS_DEM='/SNOWDATA/IDALS/REF_DEM/MCS_REFDEM_WGS84.tif'
+DC_DEM='/SNOWDATA/IDALS/REF_DEM/DC_REFDEM_WGS84.tif'
+
+ICEROAD='/SNOWDATA/IDALS/ice-road-copters/scripts/ice-road-pipeline.py'
+
+D2021='/SNOWDATA/IDALS/2021'
+D2022='/SNOWDATA/IDALS/2022'
+# SCRATCH='/home/zacharykeskinen/scratch'
+echo "Starting 2021 Flights..."
+
+for FP in $(ls $D2021)
+do
+    if [[ $FP == *"MCS"* ]]; then
+        MCS_CMD="python $ICEROAD $D2021/$FP -e $MCS_DEM -s $MCS_SHP -a $ASP_DIR"
+        echo $MCS_CMD
+        $($MCS_CMD)
+    fi
+
+    if [[ $FP == *"DC"* ]]; then
+        DC_CMD="python $ICEROAD $D2021/$FP -e $DC_DEM -s $DC_SHP -a $ASP_DIR"
+        echo $DC_CMD
+        $($DC_CMD)
+    fi
+    echo $?
+done
+
+echo $'\nStarting 2022 Flights...'
+
+for FP in $(ls $D2022)
+do
+    if [[ $FP == *"MCS"* ]]; then
+        MCS_CMD="python $ICEROAD $D2022/$FP -e $MCS_DEM -s $MCS_SHP -a $ASP_DIR"
+        echo $MCS_CMD
+        $($MCS_CMD)
+    fi
+
+    if [[ $FP == *"DC"* ]]; then
+        DC_CMD="python $ICEROAD $D2022/$FP -e $DC_DEM -s $DC_SHP -a $ASP_DIR"
+        echo $DC_CMD
+        $($DC_CMD)
+    fi
+    echo $?
+done
diff --git a/scripts/full-ice-road-run.sh b/scripts/full-ice-road-run.sh
@@ -0,0 +1,54 @@
+#!/bin/bash
+
+ASP_DIR='/SNOWDATA/IDALS/ASP/bin/'
+MCS_SHP='/SNOWDATA/IDALS/ice-road-copters/transform_area/hwy_21/hwy_21_utm_edit_v2.shp'
+DC_SHP='/SNOWDATA/IDALS/ice-road-copters/transform_areaBogusBasin/Roads_DryCreek.shp'
+MCS_DEM='/SNOWDATA/IDALS/REF_DEM/MCS_REFDEM_WGS84.tif'
+DC_DEM='/SNOWDATA/IDALS/REF_DEM/DC_REFDEM_WGS84.tif'
+
+ICEROAD='/SNOWDATA/IDALS/ice-road-copters/scripts/ice-road-pipeline.py'
+
+D2021='/SNOWDATA/IDALS/2021'
+D2022='/SNOWDATA/IDALS/2022'
+SCRATCH='/home/zacharykeskinen/scratch'
+echo "Starting 2021 Flights..."
+
+for FP in $(ls $D2021)
+do
+    $(cp -r $D2021/$FP $SCRATCH)
+    if [[ $FP == *"MCS"* ]]; then
+        MCS_CMD="python $ICEROAD $SCRATCH/$FP -e $MCS_DEM -s $MCS_SHP -a $ASP_DIR"
+        echo $MCS_CMD
+        $($MCS_CMD)
+    fi
+
+    if [[ $FP == *"DC"* ]]; then
+        DC_CMD="python $ICEROAD $SCRATCH/$FP -e $DC_DEM -s $DC_SHP -a $ASP_DIR"
+        echo $DC_CMD
+        $($DC_CMD)
+    fi
+    echo $?
+    $(cp -r $SCRATCH/$FP/ice-road $D2021/$FP)
+    $(rm -r $SCRATCH/$FP)
+done
+
+echo $'\nStarting 2022 Flights...'
+
+for FP in $(ls $D2022)
+do
+    $(cp -r $D2022/$FP $SCRATCH)
+    if [[ $FP == *"MCS"* ]]; then
+        MCS_CMD="python $ICEROAD $SCRATCH/$FP -e $MCS_DEM -s $MCS_SHP -a $ASP_DIR"
+        echo $MCS_CMD
+        $($MCS_CMD)
+    fi
+
+    if [[ $FP == *"DC"* ]]; then
+        DC_CMD="python $ICEROAD $SCRATCH/$FP -e $DC_DEM -s $DC_SHP -a $ASP_DIR"
+        echo $DC_CMD
+        $($DC_CMD)
+    fi
+    echo $?
+    $(cp -r $SCRATCH/$FP/ice-road $D2022/$FP/)
+    $(rm -r $SCRATCH/$FP)
+done