Merge pull request #82 from satsin06/extract_training

Extract LiDAR data and tests

Author: MarconiS

.github/workflows/python-package.yml

@@ -48,7 +48,7 @@ jobs:
 
     - name: Upload coverage to Codecov
       uses: codecov/codecov-action@v1
-      if: matrix.os == 'ubuntu-20.04' && matrix.python-version == '3.8' && matrix.r-version == 'release'
+      if: matrix.os == 'ubuntu-20.04' && matrix.python-version == '3.8'
       env:
         OS: ${{ runner.os }}
         PYTHON: ${{ matrix.python-version }}

neonwranglerpy/lib/extract_lidar_data.py

@@ -0,0 +1,79 @@
+"""Function to extract lidar data using RGB data and vst data."""
+import laspy
+import numpy as np
+import os
+import re
+from neonwranglerpy.lib.retrieve_aop_data import retrieve_aop_data
+
+
+def extract_lidar_data(rgb_data,
+                       vst_data,
+                       year,
+                       savepath="/content",
+                       dpID="DP1.30003.001",
+                       site="DELA"):
+    """
+    Extract LiDAR data using RGB data and vst data.
+
+    Arguments:
+    ------------
+        rgb_data: GeoDataFrame containing the plot data
+        vst_data: DataFrame containing the plot data
+        year: Year of the data
+        savepath: Path to save the data
+        dpID: LiDAR data product ID
+        site: Site name
+    """
+    retrieve_aop_data(vst_data, year, dpID, savepath)
+
+    site_level_data = vst_data[vst_data.plotID.str.contains(site)]
+    get_tiles = (((site_level_data.easting / 1000).astype(int) * 1000).astype(str) + "_" +
+                 ((site_level_data.northing / 1000).astype(int) * 1000).astype(str))
+
+    pattern = fr"{get_tiles.unique()}_classified_point_cloud.laz"
+
+    saveFile = savepath + "/" + dpID
+
+    filtered_data_list = []
+
+    for root, dirs, files in os.walk(saveFile):
+        for file in files:
+            if re.search(pattern, file):
+                lidar_file = os.path.join(root, file)
+                # directory_path = os.path.dirname(lidar_file) + '/'
+                file_name = os.path.basename(lidar_file)
+                print(lidar_file)
+
+                lidar = laspy.read(lidar_file)
+
+                x = lidar.x
+                y = lidar.y
+                z = lidar.z
+
+                data = np.vstack((x, y, z)).transpose()
+
+                lidar_dir = savepath + "/data/lidar"
+                os.makedirs(lidar_dir, exist_ok=True)
+
+                for index, row in rgb_data.iterrows():
+                    geometry = row['geometry']
+
+                    minx, miny, maxx, maxy = geometry.bounds
+
+                    filtered_data = data[(data[:, 0] >= minx) & (data[:, 0] <= maxx) &
+                                         (data[:, 1] >= miny) & (data[:, 1] <= maxy)]
+
+                    if len(filtered_data) > 0:
+                        filtered_data_list.append(filtered_data)
+
+                        filename = os.path.join(lidar_dir,
+                                                f"lidar_{file_name}_{index}.npy")
+                        np.save(filename, filtered_data)
+
+                        print(f"LiDAR data for index {index} saved as '{filename}'")
+                    else:
+                        print(f"No LiDAR data for index {index}")
+
+    filtered_data_array = np.concatenate(filtered_data_list, axis=0)
+
+    return filtered_data_array

neonwranglerpy/lib/predict_aop_data.py

@@ -4,6 +4,7 @@
 import os
 import re
 import matplotlib.pyplot as plt
+import pandas as pd
 from shapely.geometry import Point
 import rasterio
 from deepforest import main
@@ -42,8 +43,8 @@ def predict_aop_data(vst_data,
         vst_data['utmZone'].map(lambda x: (326 * 100) + int(x[:-1]))).astype(str)
     geo_data_frame = gpd.GeoDataFrame(vst_data, geometry=geometry, crs=epsg_codes.iloc[0])
     site_level_data = vst_data[vst_data.plotID.str.contains(site)]
-    get_tiles = ((site_level_data.easting / 1000).astype(int) * 1000).astype(str) + "_"
-    +((site_level_data.northing / 1000).astype(int) * 1000).astype(str)
+    get_tiles = (((site_level_data.easting / 1000).astype(int) * 1000).astype(str) + "_" +
+                 ((site_level_data.northing / 1000).astype(int) * 1000).astype(str))
     print(get_tiles.unique())
 
     pattern = fr"{year}_{site}_.*_{get_tiles.unique()[0]}"
@@ -96,4 +97,6 @@ def predict_aop_data(vst_data,
 
                         all_predictions.append(prediction)
 
-    return all_predictions
+    all_predictions_df = pd.concat(all_predictions)
+
+    return all_predictions_df

requirements.txt

@@ -14,3 +14,4 @@ opencv-python
 numpy
 rasterio
 deepforest
+laspy[lazrs,laszip]

tests/raw_data/dataframe.cpg

@@ -0,0 +1 @@
+ISO-8859-1

tests/raw_data/dataframe.dbf

@@ -0,0 +1 @@
+PROJCS["WGS_1984_UTM_Zone_16N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-87.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]

tests/raw_data/dataframe.prj

@@ -0,0 +1,22 @@
+"""Test extract_lidar_data.py file."""
+import geopandas as gpd
+import pandas as pd
+from neonwranglerpy.lib.extract_lidar_data import extract_lidar_data
+
+
+def test_extract_lidar_data():
+    """Test extract_lidar_data function."""
+    savepath = 'tests/raw_data'
+    vst_data = pd.read_csv('tests/raw_data/vst_data.csv')
+
+    rgb_data = gpd.read_file("tests/raw_data/dataframe.shp")
+
+    result = extract_lidar_data(rgb_data=rgb_data,
+                                vst_data=vst_data,
+                                year="2018",
+                                savepath=savepath,
+                                dpID="DP1.30003.001",
+                                site="DELA")
+
+    assert (vst_data.shape[0] > 0) & (vst_data.shape[1] > 0)
+    assert len(result) > 0