<jupyter_start><jupyter_text>Try to follow the tutorial at point-cloud-processing but using my own las data from NZ linz website<jupyter_code>import glob
import io
import ipyleaflet
import IPython.display
import ipyvolume.pylab as p3
import json
import matplotlib.cm
import matplotlib.pyplot as plt
import numpy as np
import os
import pandas as pd
import pdal
import PIL
import pyproj
import requests
import fiona
import shapely.geometry
import scipy.spatial
import sys
import urllib.request

%load_ext autoreload
%autoreload 2
    
sys.path.append('../src')
from pcl_utils import local_max

# Url for aerial imagery
# IVaerial = "https://geoservices.informatievlaanderen.be/raadpleegdiensten/ogw/wms?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&CRS=EPSG:31370&BBOX={0},{1},{2},{3}&WIDTH=512&HEIGHT=512&LAYERS=OGWRGB13_15VL&STYLES=default&FORMAT=image/png"

%matplotlib inline
if not os.path.isdir("../data"):
    os.makedirs("../data")
lidar_filepath = "/media/thistle/Passport/gis/nz/lidar/points.las"
lidar_filename = lidar_filepath.split("/")[-1]
if not os.path.isfile("../data/{}".format(lidar_filename)):
    subprocess.run(["cp", lidar_filepath, "../data/" + lidar_filename])
file1 = os.path.join("../data/", lidar_filename)
file1
!pdal info $file1
# arial imagery aquired from 
# https://data.linz.govt.nz/layer/51910-otago-075m-rural-aerial-photos-2004-2011/
# @ resolution 0.5
#<jupyter_output><empty_output><jupyter_text>The Data comprises:• 3,886 x ortho-rectified RGB GeoTIFF images in NZTM projection, tiled into the LINZ Standard 1:10,000 tile layout•Tile layout in NZTM projection containing relevant information.The supplied imagery is in terms of New Zealand Transverse Mercator (NZTM) map projection EPSG:2193 || NZGD2000  this dont work in Jupyter labm = ipyleaflet.Map(center=(-45.033, 168.892), zoom=12)dc = ipyleaflet.DrawControl()m.add_control(dc)m<jupyter_code># skipping the polygon clipping b/c I only have 6 million points anyway<jupyter_output><empty_output><jupyter_text>dc.last_draw  CONVERT THE LONG/LAT VALUES INTO EPSG:2193wsg84 = pyproj.Proj(init='epsg:4326')nzgd = pyproj.Proj(init="epsg:2193") coords = [pyproj.transform(wsg84, nzgd, x, y) for (x,y) in dc.last_draw['geometry']['coordinates'][0]]poly = shapely.geometry.Polygon(coords) dc.last_draw['geometry']['coordinates'][0] coords poly.wkt poly.bounds poly<jupyter_code>file1
## make a hillshade of my points.las file
file_in = "../data/points.las"
file_out = "../data/points_ground.las"
pipeline_1 = {
  "pipeline": [
      {   "type": "readers.las",
          "filename": file_in
      },
      {    "type":"filters.outlier",
           "method":"statistical",
           "mean_k":"8",
           "multiplier":"3"
      },
      {
          "type":"filters.smrf"
      },
      {
          "type":"filters.range",
          "limits":"Classification[2:2]"
      },
      {
          "type":"filters.delaunay"
      },
      {
          "type": "writers.las",
          "filename": file_out
      }
  ]}
pipeline = pdal.Pipeline(json.dumps(pipeline_1))
pipeline.validate()
%%time
np = pipeline.execute()
print(f"{np} points processed")
# write to tif and use gdal for hs
pipe2 = {
    "pipeline":[
        {
            "type":"readers.las",
            "filename":file_in
        },
        {
            "type":"writers.gdal",
            "filename":"../data/nz1.tif",
            "output_type":"all",
            "resolution":"2.0"
        }
        
    ]
}
pip3 = pdal.Pipeline(json.dumps(pipe2))
pip3.validate()
%%time
np = p3.execute()
print(f"{np} points processed")
!gdaldem hillshade ../data/nz1.tif ../data/nz1_hs.tif -z 1.0 -s 1.0 -az 315.0 -alt 45.0 -of GTiff
# !Damn looks good
# use qgis to select a polygon of interest to filter the .las file with
# maybe in this shapefile from qgis is my ticket
foi = "/home/thistle/Desktop/play.shp"
with fiona.open(foi) as f:
    temp = f.next()
    print(temp)
coords_aoi1 = temp['geometry']['coordinates'][0]
coords_aoi1
#!pdal info $file_in
p1 = shapely.geometry.Polygon(coords_aoi1)
print(p1.bounds)
# so clip 6million pt .las file by the polygon-->p1
# first try with bounds 
# next with a Polygon.wkt

t = p1.bounds
t<jupyter_output><empty_output><jupyter_text>bounds needs to be in following format"bounds":"([0,1000000],[0,1000000])"<jupyter_code>proper_format_bounds = ([t[0], t[2]], [t[1], t[3]])
str(proper_format_bounds)
!ls -al ../data
proper_format_bounds
file_in = "/home/thistle/Documents/dev/ML/gis_project/data/points.las"

file_out1 = "../data/clipped_las.las"
file_out2 = "../data/clipped_dem.tif"
file_out2
## this works for cropping by rect bounds
# see below if by polygon

pipeline_final = {
  "pipeline": [
      {   "type": "readers.las",
          "filename": file_in
      },
      {
          "type":"filters.crop",
          "bounds": str(proper_format_bounds)
      },
      {
          "type": "writers.las",
          "filename": file_out1
      },
      {    "type":"filters.outlier",
           "method":"statistical",
           "mean_k":"8",
           "multiplier":"3"
      },
      {
          "type":"filters.smrf"
      },
      {
          "type":"filters.range",
          "limits":"Classification[2:2]"
      },
      {
          "type":"filters.delaunay"
      },
      {
            "type": "writers.gdal",
            "gdaldriver": "GTIFF",
            "output_type": "all",
            "resolution": "2.0",
            "filename": file_out2
      }
  ]}
pipeline3 = pdal.Pipeline(json.dumps(pipeline_final))
pipeline3.validate()
## full pipeline
%%time
np = pipeline3.execute()
print(f"{np} pts processed")
## sweet

!gdaldem hillshade ../data/clipped_dem.tif ../data/clipped_hs.tif -z 1.0 -s 1.0 -az 315.0 -alt 45.0 -of GTiff
with fiona.open(foi) as f:
    print(f[0])
with fiona.open(foi) as f:
    temp = f[0]['geometry']['coordinates'][0]
temp
## actually -- The clipping polygon, expressed in a well-known text string, is needed
p1 = shapely.geometry.Polygon(temp)
p1.to_wkt()
## this one for polygon shape
pipeline_final_poly = {
  "pipeline": [
      {   "type": "readers.las",
          "filename": file_in
      },
      {
          "type":"filters.crop",
          "polygon": p1.to_wkt()
      },
      {
          "type": "writers.las",
          "filename": file_out1
      },
      {    "type":"filters.outlier",
           "method":"statistical",
           "mean_k":"8",
           "multiplier":"3"
      },
      {
          "type":"filters.smrf"
      },
      {
          "type":"filters.range",
          "limits":"Classification[2:2]"
      },
      {
          "type":"filters.delaunay"
      },
      {
            "type": "writers.gdal",
            "gdaldriver": "GTIFF",
            "output_type": "all",
            "resolution": "2.0",
            "filename": file_out2
      }
  ]}

pipeline_cp = pdal.Pipeline(json.dumps(pipeline_final_poly))
pipeline_cp.validate()
%%time
np = pipeline_cp.execute()
print(f"{np} points processed")
p1.to_wkt()
arr = pipeline3.arrays[0]
description = arr.dtype.descr
description
cols = [col for col, _ in description]
df=pd.DataFrame({col: arr[col] for col in cols})
arr.shape, df.shape
df.head()
# what is percentage of ground points
# (df.Classification == 2)
df.shape[0] # .78
df["X_o"] = df['X']
df['Y_o'] = df["Y"]
df["Z_o"] = df["Z"]

df['X'] = df['X'] - df['X_o'].min()
df['Y'] = df['Y'] - df['Y_o'].min()
df['Z'] = df['Z'] - df['Z_o'].min()
df.columns
t = arr[0]
print(t)
fig = p3.figure(width=1000)
fig.xlabel='Y'
fig.ylabel='Z'
fig.zlabel='X'
all_points = p3.scatter(df['Y'], df['Z'], df['X'], color='red', size=.2)
p3.squarelim()
p3.show()
# now classify grounds
file_in = ""
file_out
# practice dem by following tutorial at
# https://pdal.io/workshop/exercises/analysis/dtm/dtm.html
# first identify ground
foi = "/home/thistle/Documents/dev/ML/gis/lidar_nbs/pdal_workshop/data/exercises/analysis/ground/CSite1_orig-utm.laz"
p1 = {
  "pipeline": [
      {   "type": "readers.las",
          "filename": foi
      },
      {
          "type":"filters.smrf"
      },
      {
          "type":"filters.range",
          "limits":"Classification[2:2]"
      },
      {
          "type": "writers.las",
          "filename": "../data/ground_only.laz"
      }
  ]}
pipe1 = pdal.Pipeline(json.dumps(p1))
pipe1.validate()
%%time
np = pipe1.execute()
print(f"pipeline- pipe1 classified {np} points")<jupyter_output>pipeline- pipe1 classified 386067 points
CPU times: user 6.89 s, sys: 8.01 ms, total: 6.9 s
Wall time: 6.88 s<jupyter_text>The primary input for Digital Terrain Model generation is a point cloud with ground vs. not-ground classifications. In this example, we will use an algorithm provided by PDAL, the Simple Morphological Filter technique to generate a ground surface. In this invocation, we have more control over the process. First the outlier filter merely classifies outliers with a Classification value of 7. These outliers are then ignored during SMRF processing with the ignore option. Finally, we add a range filter to extract only the ground returns (i.e., Classification value of 2).<jupyter_code>p2 = {
  "pipeline": [
      {   "type": "readers.las",
          "filename": foi
      },
      {    "type":"filters.outlier",
           "method":"statistical",
           "mean_k":"8",
           "multiplier":"3"
      },
      {
          "type":"filters.smrf"
      },
      {
          "type":"filters.range",
          "limits":"Classification[2:2]"
      },
      {
          "type": "writers.las",
          "filename": "../data/ground_only.laz"
      }
  ]}
pipe2 = pdal.Pipeline(json.dumps(p2))
pipe2.validate()
%%time
pipe2.execute()
# with ground classified as --
infile = "../data/ground_only.laz"
outfile = "../data/dem1.tif"

p3 = {
    "pipeline":[
        {   "type":"readers.las",
            "filename": infile
        },
        {
            "type": "writers.gdal",
            "gdaldriver": "GTIFF",
            "output_type": "all",
            "resolution": "2.0",
            "filename": outfile
        }
    ]
}
pipe3 = pdal.Pipeline(json.dumps(p3))
pipe3.validate()
np = pipe3.execute()
print(f"pipe 3 processed {np} points")
p4 = {
  "pipeline": [
      {   "type": "readers.las",
          "filename": foi
      },
      {    "type":"filters.outlier",
           "method":"statistical",
           "mean_k":"8",
           "multiplier":"3"
      },
      {
          "type":"filters.smrf"
      },
      {
          "type":"filters.range",
          "limits":"Classification[2:2]"
      },
      {
          "type": "filters.poisson",
          "depth": "16"
      },
      {
          "type": "writers.las",
          "filename": "../data/ground_mesh.laz"
      }
  ]}
pipe4 = pdal.Pipeline(json.dumps(p4))
pipe4.validate()
%%time
np = pipe4.execute()
print(f"pipe4 processed {np} points")
# with ground classified as --
infile = "../data/ground_mesh.laz"
outfile = "../data/dem2_mesh.tif"

p5 = {
    "pipeline":[
        {   "type":"readers.las",
            "filename": infile
        },
        {
            "type": "writers.gdal",
            "gdaldriver": "GTIFF",
            "output_type": "all",
            "resolution": "2.0",
            "filename": outfile
        }
    ]
}
pipe5 = pdal.Pipeline(json.dumps(p5))
pipe5.validate()
np = pipe5.execute()
print(f"{np} points processed")
pipeline2 = {
  "pipeline": [
    {   "type": "readers.las",
        "filename": foi
    },
    {
        "filename":"./dtm.tif",
        "gdaldriver":"GTiff",
        "output_type":"all",
        "resolution":"2.0",
        "type": "writers.gdal"
    }
  ]
}
json.dumps(pipeline2)
pipeline2
!wget https://coast.noaa.gov/htdata/lidar1_z/geoid12b/data/5008/20131013_usgs_olympic_wa_10TDT310990.laz
!pdal info "20131013_usgs_olympic_wa_10TDT310990.laz"
import subprocess
import json
result = subprocess.call(['pdal', "info", "20131013_usgs_olympic_wa_10TDT310990.laz"])

# result.stdout.decode()
json.loads(result.stdout.decode())
!pwd
pipeline_string = """{
  "pipeline": [
    {
      "type": "readers.las",
      "filename": "/media/thistle/Passport/gis/nz/lidar/points.las"
    },
    {
      "type": "writers.text",
      "order": "Y,X,Z",
      "keep_unspecified": "false",
      "quote_header": "false",
      "delimiter": ",",
      "filename": "output1.txt"
    }
  ]
}"""

pipeline_string2 = """{
  "pipeline": [
    {
      "type": "readers.text",
      "filename": "output1.txt",
      "header": "X, Y, Z",
      "skip":1
    },
    {
      "type": "writers.text",
      "filename": "output2.txt"
      
    }
  ]
}"""
with open("../data/las2yxz.json", "w") as file:
    file.write(pipeline_string)
with open("../data/yxz2xyz.json", "w") as file2:
    file2.write(pipeline_string2)
!ls -al ../data/
# !more ../data/las2yxz.json
!pdal pipeline ../data/las2yxz.json
!pdal pipeline ../data/yxz2xyz.json
!pdal info output1.las
!cat ../data/yxz2las.json
!pdal info output1.txt<jupyter_output>{
  "filename": "output1.txt",
  "pdal_version": "2.0.1 (git-version: 6600e6)",
  "stats":
  {
    "bbox":
    {
      "native":
      {
        "bbox":
        {
          "maxx": 1277038.29,
          "maxy": 5006227.75,
          "maxz": 1404.07,
          "minx": 1275745.46,
          "miny": 5003951.5,
          "minz": 743.34
        },
        "boundary": "{ \"type\": \"Polygon\", \"coordinates\": [ [ [ 1275745.46, 5003951.5, 743.34 ], [ 1275745.46, 5006227.75, 743.34 ], [ 1277038.29, 5006227.75, 1404.07 ], [ 1277038.29, 5003951.5, 1404.07 ], [ 1275745.46, 5003951.5, 743.34 ] ] ] }"
      }
    },
    "statistic":
    [
      {
        "average": 1276561.587,
        "count": 6727923,
        "maximum": 1277038.29,
        "minimum": 1275745.46,
        "name": "X",
        "position": 0,
        "stddev": 575.6314692,
        "variance": 331351.5884
      },
      {
        "average": 5005217.851,
        "count": 6727923,
        "maximum": 5006227.75,
        "minimum": 50039[...]
{'fraction_non_alphanumeric': 0.16707404744787924, 'fraction_numerical': 0.04370956146657081, 'mean_word_length': 2.602019678922838, 'pattern_counts': {'":': 163, '<': 19, '<?xml version=': 0, '>': 20, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '556388', 'license': 'MIT', 'n_tokens_mistral': 5390, 'n_tokens_neox': 4985, 'n_words': 988, 'path': 'notebooks/nz_lidar_first_try.ipynb', 'repo_name': 'buffbob/NewZealand_LightningStrikes', 'size': '13910'}