Merge pull request #5 from ProgrammingOperative/visualization

Visualization complete

Author: ProgrammingOperative

assets/usgs_data_pipeline.json

@@ -16,7 +16,7 @@
         "type": "filters.range",
         "inputs": [ "crop" ],
         "tag": "no_noise",
-        "limits": "Classification![2:2]"
+        "limits": "Classification![7:7]"
       },
       {
         "type": "filters.assign",
@@ -25,9 +25,16 @@
         "assignment": "Classification[:]=0"
       },
       {
-        "type": "filters.smrf",
+        "type": "filters.reprojection",
         "inputs": [ "wipe_classes" ],
-        "tag": "groundify"
+        "tag": "reprojection",
+        "in_srs": "EPSG:3857",
+        "out_srs": ""
+      },
+      {
+        "inputs": [ "reprojection" ],
+        "tag": "groundify",
+        "type": "filters.smrf"
       },
       {
         "type": "filters.range",
@@ -36,11 +43,10 @@
         "limits": "Classification[2:2]"
       },
       {
-        "type": "filters.reprojection",
+        "type": "writers.las",
         "inputs": [ "classify" ],
-        "tag": "reprojection",
-        "in_srs": "EPSG:3857",
-        "out_srs": ""
+        "tag": "writers_las",
+        "filename": ""
       },
       {
         "type": "writers.gdal",

notebooks/extract_and_visualize.ipynb

@@ -4,14 +4,15 @@
 
 
 
-
-
 class AgriTechLidar:
-    def __init__(self) -> None:
-        pass
+    def __init__(self, epsg = 26915):
+        self.input_epsg = 3857
+        self.output_epsg = epsg
+        self.fetch_data = FetchLidar(self.output_epsg)
+
 
-    def fetch_lidar(self, polygon: Polygon):
-        return self.fetch_lidar(polygon)
+    def fetch_lidar(self, polygon: Polygon, regions=[]):
+        return self.fetch_data.get_lidar_data(polygon, regions)
 
 
     def render_vis(self, df: gpd.GeoDataFrame) -> Vis:

notebooks/lidar1.ipynb

@@ -5,11 +5,15 @@
 import geopandas as gpd
 from shapely.geometry import Polygon
 import json
+from gpd_helper import GpdHelper
 from config import Config
 
 class FetchLidar:
-    def __init__(self) -> None:
-        pass
+    def __init__(self, epsg: int = 26915) -> None:
+        self._input_epsg = 3857
+        self.output_epsg = epsg
+        self._gdf_helper = GPDHelper(self._input_epsg, self.output_epsg)
+
 
 
 
@@ -22,19 +26,26 @@ def make_pipeline(self, bounds: str, polygon_str: str, region: str, filename: st
         pipe['pipeline'][0]['filename'] = Config.USGS_3DEP_PUBLIC_DATA_PATH / f"{region}/etp.json"
         pipe['pipeline'][0]['bounds'] = bounds
         pipe['pipeline'][1]['polygon'] = polygon_str
-        pipe['pipeline'][6]['out_srs'] = f'EPSG:{self.output_epsg}'
-        pipe['pipeline'][7]['filename'] = Config.DATA_PATH + f'{filename}.tif'
-        pipe['pipeline'][8]['filename'] = Config.DATA_PATH + f'{filename}.tif' 
+        pipe['pipeline'][5]['out_srs'] = f'EPSG:{self.output_epsg}'
+        pipe['pipeline'][8]['filename'] = Config.DATA_PATH + f'{filename}.laz'
+        pipe['pipeline'][9]['filename'] = Config.DATA_PATH + f'{filename}.tif' 
         return pdal.Pipeline(json.dumps(pipe))
 
 
     #   def get_dep(self, bounds: Bounds, polygon_str: str, region: list) -> gpd.GeoDataFrame:
     #   We only need the polygon, not necesarily the boundaries
 
 
-
     def fetch_data(self, bounds: Bounds, polygon_str: str, region: list) -> gpd.GeoDataFrame:
-        """ Executes pdal pipeline and fetches point cloud data from a public repository.
+        filename = region + "_" + bounds.get_bound_name()
+        pl = self.make_pipeline(bounds.get_bound_str(),
+                            polygon_str, region, filename)
+        pl.execute()
+        gpd_data = self._gdf_helper.create_gdf(pl.arrays)
+        return gpd_data
+
+
+         """ Executes pdal pipeline and fetches point cloud data from a public repository.
             Using GDfHelper class creates Geopandas data frame containing geometry and elevation of the point cloud data.
 
         Args:
@@ -45,12 +56,10 @@ def fetch_data(self, bounds: Bounds, polygon_str: str, region: list) -> gpd.GeoD
         Returns:
             gpd.GeoDataFrame: Geopandas data frame containing geometry and elevation
             """
-        filename = region + "_" + bounds.get_bound_name()
-        pl = self.make_pipeline(bounds.get_bound_str(),
-                            polygon_str, region, filename)
-        pl.execute()
-        dep_data = self._gdf_helper.create_gdf(pl.arrays)
-        return dep_data
+
+    def get_lidar_data(self, polygon, regions):
+        bound, polygon_str = self._gdf_helper.get_bound_from_polygon(polygon)
+        return self.fetch_data(bound, polygon_str, regions)
 
 
 

scripts/agritechlidar.py

@@ -54,3 +54,4 @@ def create_gdf(self, array_data: np.ndarray) -> gpd.GeoDataFrame:
 
 
 
+

scripts/fetch_lidar.py

@@ -0,0 +1,60 @@
+import os
+from glob import glob
+import matplotlib.pyplot as plt
+import rasterio as rio
+from rasterio.plot import plotting_extent, show
+import geopandas as gpd
+import numpy as np
+import geopandas as gpd
+import matplotlib.image as mpimg
+import rasterio
+from rasterio.plot import show_hist
+
+
+
+class Visualize():
+    def __init__(self, df: gpd.GeoDataFrame):
+        self.df = df
+
+    
+    def plot_raster(path_to_raster):
+        """
+        displays a raster from a .tif raster file
+        args:
+            path_to_raster (str): path to the raster file
+        returns:
+            rasterio image
+        """
+        src = rasterio.open(path_to_raster)
+        fig, (axrgb, axhist) = plt.subplots(1, 2, figsize=(14,7))
+        show((src), cmap='Greys_r', contour=True, ax=axrgb)
+        show_hist(src, bins=50, histtype='stepfilled',
+                lw=0.0, stacked=False, alpha=0.3, ax=axhist)
+        plt.show()
+
+    def get_points(self):
+        """ Generates a NumPy array from point clouds data.
+        """
+        x = self.df.geometry.x
+        y = self.df.geometry.y
+        z = self.df.elevation
+        return np.vstack((x, y, z)).transpose()
+
+    def plot_2d_heatmap(df,column,title):
+        """
+        plot a 2d heat map of the terrain
+        args:
+            df (geopndas df): a geopandas dataframe demonstrating the data
+            column (str): input column to outline in string
+            title (str): input title of the map in string
+        return:
+            2d heat map of terrain
+        """
+        fig, ax = plt.subplots(1, 1, figsize=(12, 10))
+        fig.patch.set_alpha(0)
+        plt.grid('on', zorder=0)
+        df.plot(column=column, ax=ax, legend=True, cmap="terrain")
+        plt.title(title)
+        plt.xlabel('long')
+        plt.ylabel('lat')
+        plt.show()

scripts/gpd_helper.py

@@ -0,0 +1,39 @@
+import unittest
+import sys, os
+sys.path.append(os.path.abspath(os.path.join('..')))
+
+
+class TestCases(unittest.TestCase):
+    def test_find_average(self):
+        """
+        Test that it retunrs the average of a given list
+        """
+        data = [1, 2, 3]
+        result = loader.find_average(data)
+        self.assertEqual(result, 2.0)
+
+    def test_input_value(self):
+        """
+        Provide an assertion level for arg input
+        """
+        
+        self.assertRaises(TypeError, loader.find_average, True)
+
+class TestCountOccurence(unittest.TestCase):
+    def test_count_occurence(self):
+        """
+        Test that it returns the count of each unique values in the given list
+        """
+        data = [0,0,9,0,8,9,0,7]
+        result = loader.count_occurence(data)
+        output = {0: 4, 9: 2, 8: 1, 7: 1}
+        self.assertAlmostEqual(result, output)
+
+    def test_input_value(self):
+        """
+        Provide an assertion level for arg input
+        """
+        self.assertRaises(TypeError, loader.count_occurence, True)
+
+if __name__ == '__main__':
+    unittest.main()