<jupyter_start><jupyter_text>LiDAR point cloud of the Trail Islands in British Columbia, CanadaThis is a LiDAR point cloud from the [LidarBC](https://www2.gov.bc.ca/gov/content/data/geographic-data-services/lidarbc) project, sliced to the small [Trail Islands](https://apps.gov.bc.ca/pub/bcgnws/names/21973.html) North of Vancouver to reduce the data size. The islands have some nice looking topography and their isolated nature creates problems for some interpolation methods.Coordinates of the original data are in NAD83 UTM zone 10 and the vertical datum for elevation is CGVD2013 (based on the Canadian geoid).The horizontal datum was changed to WGS84 and coordinates un-projected to longitude, latitude.License: [Open Government Licence - British Columbia](https://www2.gov.bc.ca/gov/content/data/open-data/open-government-licence-bc)<jupyter_code>import os
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import verde as vd
import pooch
import pyproj
import pygmt
import laspy<jupyter_output><empty_output><jupyter_text>Download the tile we needThe islands are covered by the tile `bc_092g041` of the dataset. Download it with Pooch and cache it.<jupyter_code>fname = pooch.retrieve(
    url="https://nrs.objectstore.gov.bc.ca/gdwuts/092/092g/2019/pointcloud/bc_092g041_4_2_2_xyes_8_utm10_2019.laz",
    known_hash="md5:5acad29ff5ba1f2a09453ed16c3004da",
)
print(f"size: {os.path.getsize(fname) / 1e6} Mb")<jupyter_output>size: 83.871764 Mb<jupyter_text>Read the dataUse laspy to read in the data from the file. We'll read it in chunks and only take the ground reflectance points to avoid overloading our memory.<jupyter_code>easts, norths, heights = [], [], []
with laspy.open(fname) as lasfile:
    for chunk in lasfile.chunk_iterator(100_000):
        ground = chunk.classification == 2
        easts.append(chunk.x.scaled_array()[ground].copy())
        norths.append(chunk.y.scaled_array()[ground].copy())
        heights.append(chunk.z.scaled_array()[ground].copy())
        
data_full = pd.DataFrame({
    "easting": np.concatenate(easts),
    "northing": np.concatenate(norths),
    "height": np.concatenate(heights),
})
data_full<jupyter_output><empty_output><jupyter_text>Cut the dataSlice the data to just the islands to avoid having too large a dataset.<jupyter_code>inside = vd.inside((data_full.easting, data_full.northing), region=(440600, 441600, 5478100, 5479000))
data = data_full[inside].reset_index(drop=True)
data<jupyter_output><empty_output><jupyter_text>Plot the point cloudMake a quick plot to make sure the data look OK.<jupyter_code>plt.figure(figsize=(15, 12))
plt.scatter(data.easting, data.northing, c=data.height, s=1)
plt.colorbar(label="meters", pad=0, aspect=50)
plt.axis("scaled")
plt.xlabel("UTM Easting (m)")
plt.ylabel("UTM Northing (m)")
plt.savefig("preview.jpg", dpi=200, bbox_inches="tight")<jupyter_output><empty_output><jupyter_text>Convert UTM to latitude, longitudeUndo the projection and store the latitude, longitude so we don't need to preserve the UTM zone information.<jupyter_code>wgs84 = pyproj.CRS.from_epsg(4326)
utm_nad83 = pyproj.CRS.from_proj4("+proj=utm +zone=10 +north +ellps=GRS80 +datum=NAD83 +units=m +no_defs")
transformer = pyproj.Transformer.from_crs(utm_nad83, wgs84, always_xy=True)
data["longitude"], data["latitude"] = transformer.transform(data.easting, data.northing)
plt.figure(figsize=(15, 12))
plt.scatter(data.longitude, data.latitude, c=data.height, s=1)
plt.colorbar(label="meters", pad=0, aspect=50)
plt.xlabel("Longitude")
plt.ylabel("Latitude")<jupyter_output><empty_output><jupyter_text>ExportMake a separate DataFrame to export to CSV with centimeter level precision.<jupyter_code>export = pd.DataFrame({
    "longitude": data.longitude.map(lambda x: "{:.7f}".format(x)),
    "latitude": data.latitude.map(lambda x: "{:.7f}".format(x)),    
    "elevation_m": data.height.map(lambda x: "{:.2f}".format(x)),
})
export
fname = "british-columbia-lidar.csv.xz"
export.to_csv(fname, index=False)
print(fname)
print(f"size: {os.path.getsize(fname) / 1e6} Mb")
for alg in ["md5", "sha256"]:
    print(f"{alg}:{pooch.file_hash(fname, alg=alg)}")<jupyter_output>british-columbia-lidar.csv.xz
size: 4.378076 Mb
md5:354c725a95036bd8340bc14e043ece5a
sha256:03c6f1b99374b8c00c424c788cb6956bc00ab477244bb69835d4171312714fe1<jupyter_text>Read back the data and plot itVerify that the output didn't corrupt anything.<jupyter_code>lidar = pd.read_csv(fname)
lidar
fig = pygmt.Figure()
fig.basemap(
    region=vd.pad_region(vd.get_region((lidar.longitude, lidar.latitude)), 0.1),
    projection="M15c",
    frame=True,
)
pygmt.makecpt(cmap="viridis", series=[lidar.elevation_m.min(), lidar.elevation_m.max()])
fig.plot(x=lidar.longitude, y=lidar.latitude, color=lidar.elevation_m, cmap=True, style="c0.1c")
fig.colorbar(frame="af")
fig.coast(shorelines=True, resolution="f")
fig.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13016907720513343, 'fraction_numerical': 0.04909350173151355, 'mean_word_length': 6.28486646884273, 'pattern_counts': {'":': 6, '<': 31, '<?xml version=': 0, '>': 31, 'https://': 4, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '122867', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 1984, 'n_tokens_neox': 1791, 'n_words': 383, 'path': 'british-columbia-lidar/prepare.ipynb', 'repo_name': 'fatiando/rockhound-data', 'size': '4909'}