1
1
"""Module for PointCloud class."""
2
+
2
3
from __future__ import annotations
3
4
4
5
import os .path
5
- import warnings
6
- from typing import Iterable , Literal , Any
7
6
import pathlib
7
+ import warnings
8
+ from typing import Any , Iterable , Literal
8
9
10
+ import geopandas as gpd
11
+ import numpy as np
9
12
import pandas as pd
10
13
from pyproj import CRS
11
- import numpy as np
12
- import geopandas as gpd
13
14
from rasterio .coords import BoundingBox
14
15
from shapely .geometry .base import BaseGeometry
15
16
16
17
import geoutils as gu
18
+ from geoutils ._typing import ArrayLike , NDArrayNum , Number
17
19
from geoutils .interface .gridding import _grid_pointcloud
18
- from geoutils .raster .sampling import subsample_array
19
20
from geoutils .interface .raster_point import _raster_to_pointcloud
20
- from geoutils ._typing import Number , NDArrayNum , ArrayLike
21
+ from geoutils .raster . sampling import subsample_array
21
22
22
23
try :
23
24
import laspy
25
+
24
26
has_laspy = True
25
27
except ImportError :
26
28
has_laspy = False
@@ -36,15 +38,18 @@ def _load_laspy_data(filename: str, columns: list[str]) -> gpd.GeoDataFrame:
36
38
data = np .vstack ([las [n ] for n in columns ]).T
37
39
38
40
# Build geodataframe
39
- gdf = gpd .GeoDataFrame (geometry = gpd .points_from_xy (x = las .x , y = las .y ,
40
- crs = las .header .parse_crs (prefer_wkt = False )),
41
- data = data ,
42
- columns = columns )
41
+ gdf = gpd .GeoDataFrame (
42
+ geometry = gpd .points_from_xy (x = las .x , y = las .y , crs = las .header .parse_crs (prefer_wkt = False )),
43
+ data = data ,
44
+ columns = columns ,
45
+ )
43
46
44
47
return gdf
45
48
46
49
47
- def _load_laspy_metadata (filename : str , ) -> tuple [CRS , int , BoundingBox , pd .Index ]:
50
+ def _load_laspy_metadata (
51
+ filename : str ,
52
+ ) -> tuple [CRS , int , BoundingBox , pd .Index ]:
48
53
"""Load point cloud metadata from LAS/LAZ/COPC file."""
49
54
50
55
with laspy .open (filename ) as f :
@@ -62,7 +67,7 @@ def _load_laspy_metadata(filename: str, ) -> tuple[CRS, int, BoundingBox, pd.Ind
62
67
#
63
68
# with laspy.open(filename) as f:
64
69
# new_hdr = laspy.LasHeader(version="1.4", point_format=6)
65
- # # You can set the scales and offsets to values tha suits your data
70
+ # # You can set the scales and offsets to values that suits your data
66
71
# new_hdr.scales = np.array([1.0, 0.5, 0.1])
67
72
# new_las = laspy.LasData(header = new_hdr, points=)
68
73
#
@@ -91,9 +96,11 @@ class PointCloud(gu.Vector):
91
96
See the API for more details.
92
97
"""
93
98
94
- def __init__ (self ,
95
- filename_or_dataset : str | pathlib .Path | gpd .GeoDataFrame | gpd .GeoSeries | BaseGeometry ,
96
- data_column : str | None = "z" ):
99
+ def __init__ (
100
+ self ,
101
+ filename_or_dataset : str | pathlib .Path | gpd .GeoDataFrame | gpd .GeoSeries | BaseGeometry ,
102
+ data_column : str | None = "z" ,
103
+ ):
97
104
"""
98
105
Instantiate a point cloud from either a data column name and a vector (filename, GeoPandas dataframe or series,
99
106
or a Shapely geometry), or only with a point cloud file type.
@@ -117,8 +124,10 @@ def __init__(self,
117
124
return
118
125
# For filename, rely on parent Vector class or LAS file reader
119
126
else :
120
- if isinstance (filename_or_dataset , (str , pathlib .Path )) and \
121
- os .path .splitext (filename_or_dataset )[- 1 ] in [".las" , ".laz" ]:
127
+ if isinstance (filename_or_dataset , (str , pathlib .Path )) and os .path .splitext (filename_or_dataset )[- 1 ] in [
128
+ ".las" ,
129
+ ".laz" ,
130
+ ]:
122
131
# Load only metadata, and not the data
123
132
crs , nb_points , bounds , columns = _load_laspy_metadata (filename_or_dataset )
124
133
self ._name = filename_or_dataset
@@ -132,19 +141,18 @@ def __init__(self,
132
141
super ().__init__ (filename_or_dataset )
133
142
# Verify that the vector can be cast as a point cloud
134
143
if not all (p == "Point" for p in self .ds .geom_type ):
135
- raise ValueError ("This vector file contains non-point geometries, "
136
- "cannot be instantiated as a point cloud." )
144
+ raise ValueError (
145
+ "This vector file contains non-point geometries, " "cannot be instantiated as a point cloud."
146
+ )
137
147
138
148
# Set data column following user input
139
149
self .set_data_column (new_data_column = data_column )
140
150
141
-
142
151
# TODO: Could also move to Vector directly?
143
152
##############################################
144
153
# OVERRIDDEN VECTOR METHODS TO SUPPORT LOADING
145
154
##############################################
146
155
147
-
148
156
@property
149
157
def ds (self ) -> gpd .GeoDataFrame :
150
158
"""Geodataframe of the point cloud."""
@@ -168,15 +176,15 @@ def ds(self, new_ds: gpd.GeoDataFrame | gpd.GeoSeries) -> None:
168
176
@property
169
177
def crs (self ) -> CRS :
170
178
"""Coordinate reference system of the vector."""
171
- # Overridding method in Vector
179
+ # Overriding method in Vector
172
180
if self .is_loaded :
173
181
return super ().crs
174
182
else :
175
183
return self ._crs
176
184
177
185
@property
178
186
def bounds (self ) -> BoundingBox :
179
- # Overridding method in Vector
187
+ # Overriding method in Vector
180
188
if self .is_loaded :
181
189
return super ().bounds
182
190
else :
@@ -185,7 +193,7 @@ def bounds(self) -> BoundingBox:
185
193
@property
186
194
def all_columns (self ) -> pd .Index :
187
195
"""Index of all columns of the point cloud, excluding the column of 2D point geometries."""
188
- # Overridding method in Vector
196
+ # Overriding method in Vector
189
197
if self .is_loaded :
190
198
all_columns = super ().columns
191
199
all_columns_nongeom = all_columns [all_columns != "geometry" ]
@@ -210,8 +218,10 @@ def set_data_column(self, new_data_column: str) -> None:
210
218
"""Set new column as point cloud data column."""
211
219
212
220
if new_data_column not in self .all_columns :
213
- raise ValueError (f"Data column { new_data_column } not found among columns, available columns "
214
- f"are: { ', ' .join (self .all_columns )} ." )
221
+ raise ValueError (
222
+ f"Data column { new_data_column } not found among columns, available columns "
223
+ f"are: { ', ' .join (self .all_columns )} ."
224
+ )
215
225
self ._data_column = new_data_column
216
226
217
227
@property
@@ -266,13 +276,16 @@ def from_array(cls, array: NDArrayNum, crs: CRS, data_column: str | None = "z")
266
276
array_in = array
267
277
268
278
# Build geodataframe
269
- gdf = gpd .GeoDataFrame (geometry = gpd .points_from_xy (x = array_in [0 , :], y = array_in [1 , :], crs = crs ),
270
- data = {data_column : array_in [2 , :]})
279
+ gdf = gpd .GeoDataFrame (
280
+ geometry = gpd .points_from_xy (x = array_in [0 , :], y = array_in [1 , :], crs = crs ), data = {data_column : array_in [2 , :]}
281
+ )
271
282
272
283
return cls (filename_or_dataset = gdf , data_column = data_column )
273
284
274
285
@classmethod
275
- def from_tuples (cls , tuples_xyz : Iterable [tuple [Number , Number , Number ]], crs : CRS , data_column : str | None = "z" ) -> PointCloud :
286
+ def from_tuples (
287
+ cls , tuples_xyz : Iterable [tuple [Number , Number , Number ]], crs : CRS , data_column : str | None = "z"
288
+ ) -> PointCloud :
276
289
"""Create point cloud from an iterable of 3-tuples (X coordinate, Y coordinate, Z value)."""
277
290
278
291
return cls .from_array (np .array (tuples_xyz ), crs = crs , data_column = data_column )
@@ -316,30 +329,40 @@ def pointcloud_equal(self, other: PointCloud, **kwargs: Any):
316
329
317
330
return all ([vector_eq , data_column_eq ])
318
331
319
- def grid (self ,
320
- ref : gu .Raster | None ,
321
- grid_coords : tuple [np .ndarray , np .ndarray ] | None ,
322
- resampling : Literal ["nearest" , "linear" , "cubic" ],
323
- dist_nodata_pixel : float = 1. ) -> gu .Raster :
332
+ def grid (
333
+ self ,
334
+ ref : gu .Raster | None ,
335
+ grid_coords : tuple [np .ndarray , np .ndarray ] | None ,
336
+ resampling : Literal ["nearest" , "linear" , "cubic" ],
337
+ dist_nodata_pixel : float = 1.0 ,
338
+ ) -> gu .Raster :
324
339
"""Grid point cloud into a raster."""
325
340
326
341
if isinstance (ref , gu .Raster ):
327
342
if grid_coords is None :
328
- warnings .warn ("Both reference raster and grid coordinates were passed for gridding, "
329
- "using only the reference raster." )
343
+ warnings .warn (
344
+ "Both reference raster and grid coordinates were passed for gridding, "
345
+ "using only the reference raster."
346
+ )
330
347
grid_coords = ref .coords (grid = False )
331
348
else :
332
349
grid_coords = grid_coords
333
350
334
- array , transform = _grid_pointcloud (self .ds , grid_coords = grid_coords , data_column_name = self .data_column ,
335
- resampling = resampling , dist_nodata_pixel = dist_nodata_pixel )
351
+ array , transform = _grid_pointcloud (
352
+ self .ds ,
353
+ grid_coords = grid_coords ,
354
+ data_column_name = self .data_column ,
355
+ resampling = resampling ,
356
+ dist_nodata_pixel = dist_nodata_pixel ,
357
+ )
336
358
337
359
return gu .Raster .from_array (data = array , transform = transform , crs = self .crs , nodata = None )
338
360
339
361
def subsample (self , subsample : float | int , random_state : int | np .random .Generator | None = None ) -> PointCloud :
340
362
341
- subsample = subsample_array (array = self .ds [self .data_column ].values , subsample = subsample ,
342
- return_indices = True , random_state = random_state )
363
+ subsample = subsample_array (
364
+ array = self .ds [self .data_column ].values , subsample = subsample , return_indices = True , random_state = random_state
365
+ )
343
366
344
367
return PointCloud (self .ds [subsample ])
345
368
0 commit comments