Access CMS LiDAR data

import earthaccess
import numpy as np
import pandas as pd
import geopandas as gpd
import contextily as cx
import laspy
from shapely.geometry import box
from shapely.ops import orient

# works if the EDL login already been persisted to a netrc
auth = earthaccess.login(strategy="netrc") 
if not auth.authenticated:
    # ask for EDL credentials and persist them in a .netrc file
    auth = earthaccess.login(strategy="interactive", persist=True)

granules = earthaccess.search_data(
    count=-1, # needed to retrieve all granules
    granule_name = '*.laz', # only find laz files
    doi="10.3334/ORNLDAAC/1644" # data collection doi
)
print(f"Total granules found: {len(granules)}")

granules[0]

def convert_umm_bbox(b):
    """converts UMM bounding rectangles"""
    try:
        return box(b[0]['WestBoundingCoordinate'], 
              b[0]['SouthBoundingCoordinate'],
              b[0]['EastBoundingCoordinate'], 
              b[0]['NorthBoundingCoordinate'])
    except TypeError as e:
        return None

def convert_list_gdf(datag):
    """converts List[] to geopandas dataframe"""
    # create pandas dataframe from json
    df = pd.json_normalize([vars(granule)['render_dict'] for granule in datag])
    # keep only last string of the column names
    df.columns=df.columns.str.split('.').str[-1]
    df['geometry'] = df['BoundingRectangles'].apply(convert_umm_bbox)
    return gpd.GeoDataFrame(df, geometry="geometry", crs="EPSG:4326")
    
gdf = convert_list_gdf(granules)

gdf[['GranuleUR', 'size', 'geometry']].explore(color='red',  fillcolor='red')

poly = gpd.read_file("assets/reserva_ducke.json")
poly.explore(color='red',  fill=False)

# bounding lon, lat as a list of tuples
poly.geometry = poly.geometry.apply(orient, args=(1,))
# simplifying the polygon to bypass the coordinates 
# limit of the CMR with a tolerance of .005 degrees
xy = poly.geometry.simplify(0.005).get_coordinates()

granules = earthaccess.search_data(
    count=-1, # needed to retrieve all granules
    granule_name = '*.laz', # only find laz files
    doi="10.3334/ORNLDAAC/1644", # dataset doi
    polygon=list(zip(xy.x, xy.y))
)
print(f"Total granules found: {len(granules)}")

gdf = convert_list_gdf(granules)
gdf[['GranuleUR', 'size', 'geometry']].explore(color='red',  fillcolor='red')

downloaded_files = earthaccess.download(granules, local_path="downloads")

downloaded_files

point_cloud = laspy.read(downloaded_files[0])
print(point_cloud)

point_format = point_cloud.point_format
list(point_format.dimension_names)