From 8fd7619e83d74e948ca6206382dd2ee531832d45 Mon Sep 17 00:00:00 2001 From: thijsvl Date: Tue, 15 Aug 2023 11:03:37 +0200 Subject: [PATCH 01/11] Very first proof of concept for Morton order sorting of a numpy array. Note that for application to an stm, we should still check whether we can specify to mortonorder() the dimensions to sort on. Otherwise, we could temporarily rename lat-lon to Y-X before sorting. --- morton_order.ipynb | 2067 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 2067 insertions(+) create mode 100644 morton_order.ipynb diff --git a/morton_order.ipynb b/morton_order.ipynb new file mode 100644 index 0000000..ad071f7 --- /dev/null +++ b/morton_order.ipynb @@ -0,0 +1,2067 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "8d3a068c-9860-40b8-b7ee-806f5eea2879", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:    (points: 318225, time: 10)\n",
+       "Coordinates:\n",
+       "    azimuth    (points) int64 dask.array<chunksize=(5000,), meta=np.ndarray>\n",
+       "    lat        (points) float32 dask.array<chunksize=(5000,), meta=np.ndarray>\n",
+       "    lon        (points) float32 dask.array<chunksize=(5000,), meta=np.ndarray>\n",
+       "    range      (points) int64 dask.array<chunksize=(5000,), meta=np.ndarray>\n",
+       "  * time       (time) int64 0 1 2 3 4 5 6 7 8 9\n",
+       "Dimensions without coordinates: points\n",
+       "Data variables:\n",
+       "    amplitude  (points, time) float32 dask.array<chunksize=(5000, 10), meta=np.ndarray>\n",
+       "    complex    (points, time) complex64 dask.array<chunksize=(5000, 10), meta=np.ndarray>\n",
+       "    phase      (points, time) float32 dask.array<chunksize=(5000, 10), meta=np.ndarray>
" + ], + "text/plain": [ + "\n", + "Dimensions: (points: 318225, time: 10)\n", + "Coordinates:\n", + " azimuth (points) int64 dask.array\n", + " lat (points) float32 dask.array\n", + " lon (points) float32 dask.array\n", + " range (points) int64 dask.array\n", + " * time (time) int64 0 1 2 3 4 5 6 7 8 9\n", + "Dimensions without coordinates: points\n", + "Data variables:\n", + " amplitude (points, time) float32 dask.array\n", + " complex (points, time) complex64 dask.array\n", + " phase (points, time) float32 dask.array" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import xarray as xr\n", + "import stmtools\n", + "\n", + "stmat = xr.open_zarr('../data/stm.zarr')\n", + "stmat" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "20ed67f7-82fa-4bcc-86fb-7353c60a8aa3", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:    (points: 316762, time: 10)\n",
+       "Coordinates:\n",
+       "    azimuth    (points) int64 dask.array<chunksize=(5000,), meta=np.ndarray>\n",
+       "    lat        (points) float32 dask.array<chunksize=(5000,), meta=np.ndarray>\n",
+       "    lon        (points) float32 dask.array<chunksize=(5000,), meta=np.ndarray>\n",
+       "    range      (points) int64 dask.array<chunksize=(5000,), meta=np.ndarray>\n",
+       "  * time       (time) int64 0 1 2 3 4 5 6 7 8 9\n",
+       "Dimensions without coordinates: points\n",
+       "Data variables:\n",
+       "    amplitude  (points, time) float32 dask.array<chunksize=(5000, 10), meta=np.ndarray>\n",
+       "    complex    (points, time) complex64 dask.array<chunksize=(5000, 10), meta=np.ndarray>\n",
+       "    phase      (points, time) float32 dask.array<chunksize=(5000, 10), meta=np.ndarray>
" + ], + "text/plain": [ + "\n", + "Dimensions: (points: 316762, time: 10)\n", + "Coordinates:\n", + " azimuth (points) int64 dask.array\n", + " lat (points) float32 dask.array\n", + " lon (points) float32 dask.array\n", + " range (points) int64 dask.array\n", + " * time (time) int64 0 1 2 3 4 5 6 7 8 9\n", + "Dimensions without coordinates: points\n", + "Data variables:\n", + " amplitude (points, time) float32 dask.array\n", + " complex (points, time) complex64 dask.array\n", + " phase (points, time) float32 dask.array" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "subset = stmat.stm.subset(method='threshold', var='amplitude', threshold='>1')\n", + "subset" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "4490ac1d-9a67-497d-8dda-a9b508e9054e", + "metadata": {}, + "outputs": [], + "source": [ + "#pipeline = pdal.Filter.mortonorder().pipeline(subset)\n", + "#pipeline = pdal.Filter.sort(dimension=\"lat\").pipeline(subset.to_array())\n", + "#print(pipeline.execute())" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "59996e0d-ed18-4182-940e-4ee2cf3b530a", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "4004326\n" + ] + } + ], + "source": [ + "import pdal\n", + "\n", + "data = \"/storage/nD_PointClouds/data/entwine/red-rocks-mpi.laz\"\n", + "\n", + "pipeline_read = pdal.Reader.las(filename=data).pipeline()\n", + "print(pipeline_read.execute()) # 1065 points" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "b2567f34-cc56-4723-899c-ac1d78b6ba3b", + "metadata": {}, + "outputs": [], + "source": [ + "# Get the data from the first array\n", + "# [array([(637012.24, 849028.31, 431.66, 143, 1,\n", + "# 1, 1, 0, 1, -9., 132, 7326, 245380.78254963, 68, 77, 88),\n", + "# dtype=[('X', ' 1940\n", + "z = arr[arr[\"Z\"] > 1960]\n", + "print(len(z)) # 704 points" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "9889abdb-ef7b-4a4e-8bf5-2cf9663a0a59", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "153927\n" + ] + } + ], + "source": [ + "# Now use pdal to clamp points that have z 1970 <= v < 1980\n", + "pipeline_clamp = pdal.Filter.range(limits=\"Z[1970:1980)\").pipeline(z)\n", + "print(pipeline_clamp.execute()) # 387 points\n", + "clamped = pipeline_clamp.arrays[0]" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "b634b515-7577-4ea4-9f1e-87308b2089ad", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "4004326\n" + ] + } + ], + "source": [ + "pipeline_sort = pdal.Filter.sort(dimension=\"X\").pipeline(arr)\n", + "print(pipeline_sort.execute())" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "3402cda4-de2e-4186-84ad-c05230995f50", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "4004326\n" + ] + } + ], + "source": [ + "pipeline_morton = pdal.Filter.mortonorder().pipeline(arr)\n", + "print(pipeline_morton.execute())" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "af3eac5a-8301-4596-a1cf-29995abbf924", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 0, 1, 2, ..., 6627, 6628, 6629], dtype=uint16)" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import numpy as np\n", + "\n", + "arr_morton = pipeline_morton.arrays[0]\n", + "arr_morton[\"Intensity\"] = np.array(range(0, len(arr_morton))) % 65536\n", + "arr_morton[\"Intensity\"]" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "0a03143a-bb03-41b9-8608-d8665db67359", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "4004326\n" + ] + } + ], + "source": [ + "output_data = \"/storage/nD_PointClouds/data/entwine/red-rocks-mpi_morton.laz\"\n", + "pipeline = pdal.Writer.las(\n", + " filename=output_data,\n", + " offset_x=\"auto\",\n", + " offset_y=\"auto\",\n", + " offset_z=\"auto\",\n", + " scale_x=0.01,\n", + " scale_y=0.01,\n", + " scale_z=0.01,\n", + ").pipeline(arr_morton)\n", + "print(pipeline.execute())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4da9327b-d6b3-4ef0-b67f-8fde0f8aa2a8", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.4" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From 4b282772dcb99f744cc988bbb3d6ab3ff74be2c7 Mon Sep 17 00:00:00 2001 From: thijsvl Date: Thu, 17 Aug 2023 11:29:07 +0200 Subject: [PATCH 02/11] Initial working example on zarr data. --- morton_order.ipynb | 751 ++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 666 insertions(+), 85 deletions(-) diff --git a/morton_order.ipynb b/morton_order.ipynb index ad071f7..50d92d7 100644 --- a/morton_order.ipynb +++ b/morton_order.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": 158, "id": "8d3a068c-9860-40b8-b7ee-806f5eea2879", "metadata": {}, "outputs": [ @@ -384,7 +384,7 @@ "Data variables:\n", " amplitude (points, time) float32 dask.array<chunksize=(5000, 10), meta=np.ndarray>\n", " complex (points, time) complex64 dask.array<chunksize=(5000, 10), meta=np.ndarray>\n", - " phase (points, time) float32 dask.array<chunksize=(5000, 10), meta=np.ndarray>