Open hcwinsemius opened 2 months ago
The question is how to do this in a sufficiently performant way. The first thing I'll try is pdal's colorization filter which allows us to attribute from a gdal supported raster to any attribute in our point cloud: https://pdal.io/en/2.7-maintenance/stages/filters.colorization.html#filters-colorization
In this way, we can use grid ID's to attribute to the point cloud, and then use that attribution as our tie-in to extract statistics from the point cloud based on those IDs.
Then we can use filter.stats to extract whichever statistic we are interested in using a where clause to isolate the grid cell in question. https://pdal.io/en/2.7-maintenance/stages/filters.stats.html#filters-stats
The question is how to do this in a sufficiently performant way. The first thing I'll try is pdal's colorization filter which allows us to attribute from a gdal supported raster to any attribute in our point cloud: https://pdal.io/en/2.7-maintenance/stages/filters.colorization.html#filters-colorization
In this way, we can use grid ID's to attribute to the point cloud, and then use that attribution as our tie-in to extract statistics from the point cloud based on those IDs.
https://github.com/localdevices/ORProfile/blob/optimizer/orprofile%2Fapi%2Fdepth_sample.py does a similar resampling approach for the sonar samples. Perhaps components can be reused?
A bit of progress. I've got stuff I might want to do differently, but keeping it simple and proof of concept first.
As such, I'm just sub-sampling the dry bed point cloud first and writing to a csv, and then feeding that in. Once everything is working, I'll likely integrate this further as with PDAL's python API and/or skip writing of an intermediate file. But in the meantime, we have a json file to tell pdal what to do:
[
"input.laz",
{
"type":"filters.sample",
"radius":2
},
{
"type":"writers.text",
"format":"csv",
"order":"X,Y,Z",
"keep_unspecified":"false",
"filename":"filtered.csv"
}
]
And can invoke that as so:
pdal pipeline filter.json
Now we have a csv that we can pull in.
@hcwinsemius -- can you show me how you generate the grid for Bamboi? I need to generate that directly to fully test (rather than a geojson) so that it has all the appropriate properties. Here's my attempt in importing:
import geopandas as gpd
import pandas as pd
import numpy as np
import orprofile.api.mesh
import numpy as np
def regrid_dry(mesh, gdf, nodata=-9999, column="Depth"):
"""
regrid a set of depth samples to the provided mesh object
Parameters
----------
mesh : orprofile.Mesh
Returns
-------
"""
# get the index numbers of each mesh cell
da_index = mesh._get_index_da()
# get index number of each sample
points = da_index.ugrid.sel_points(
x=gdf.geometry.x,
y=gdf.geometry.y
)
face_index = np.ones(len(gdf), dtype=np.int64) * nodata # -9999 is the nodata value
face_index[points.mesh2d_index] = points.values
gdf["face_index"] = face_index
depth_mean = gdf[gdf["face_index"] != nodata][[column, "face_index"]].groupby("face_index").mean()
depth_grid = np.zeros(da_index.shape) * np.nan
depth_grid[depth_mean.index] = depth_mean[column]
ds = mesh._get_empty_ds()
ds["samples"] = ("mesh2d_nFaces", depth_grid)
return ds
df = pd.read_csv('filtered.csv')
gdf = gpd.GeoDataFrame(
df, geometry=gpd.points_from_xy(df.X, df.Y), crs="EPSG:32630"
)
print(gdf.head())
mesh = gpd.read_file("mesh_example_v2.geojson")
print(mesh.head())
regrid_dry(mesh, gdf)
Which errors on out __getattr__
:
python dry_sample.py
X Y Z geometry
0 606773.98 901452.56 120.49 POINT (606773.98 901452.56)
1 606775.50 901451.27 120.07 POINT (606775.5 901451.27)
2 606776.91 901449.80 119.82 POINT (606776.91 901449.8)
3 606785.67 901435.10 120.65 POINT (606785.67 901435.1)
4 606785.28 901433.41 121.65 POINT (606785.28 901433.41)
mesh2d_nFaces yi ... cols geometry
0 0 0.000000 ... 0 POLYGON ((606783.628 901576.298, 606776.657 90...
1 1 2.050305 ... 0 POLYGON ((606784.203 901574.373, 606777.264 90...
2 2 4.097750 ... 0 POLYGON ((606784.778 901572.449, 606777.872 90...
3 3 6.141944 ... 0 POLYGON ((606785.352 901570.524, 606778.48 901...
4 4 8.182641 ... 0 POLYGON ((606785.927 901568.599, 606779.088 90...
[5 rows x 9 columns]
Traceback (most recent call last):
File "/home/smathermather/Documents/OpenDroneMap/TAHMO/Bamboi/average/dry_sample.py", line 46, in <module>
regrid_samples(mesh, gdf)
File "/home/smathermather/Documents/OpenDroneMap/TAHMO/Bamboi/average/dry_sample.py", line 20, in regrid_samples
da_index = mesh._get_index_da()
File "/home/smathermather/Documents/OpenDroneMap/TAHMO/Bamboi/average/.venv/lib/python3.10/site-packages/pandas/core/generic.py", line 6299, in __getattr__
return object.__getattribute__(self, name)
AttributeError: 'GeoDataFrame' object has no attribute '_get_index_da'