The georeferenced point cloud (PointCloud)

The georeferenced point cloud (`PointCloud`)#

A point cloud represents 2D point geometries of georeferenced coordinates associated with a main 1D data array, and optionally auxiliary data.

Although a subtype of Vector, point clouds have a very different nature than other vectors and are ubiquitous in geospatial analysis, requiring their own object type. For numerical operations, a point cloud manipulation is facilitated by its own arithmetic to manipulate its main data array, as well as a specific interface with rasters (interpolation or reduction to same coordinates, or point gridding), or other point clouds (pairwise distance matching). For geometric operations, point clouds interface in specific ways with other vector-types (zonal statistics, geometric masking).

GeoUtils aims to support these features, as well the reading and writing of point clouds both from vector-type files (e.g., ESRI shapefile, geopackage, geoparquet) usually used for sparse point clouds, and from point-cloud-type files (e.g., LAS, LAZ, COPC) usually used for dense point clouds.

Warning

Support for LAS files is still preliminary and loads all data in memory for most operations. We are working on adding operations with chunked reading.

Below, a summary of the PointCloud object and its methods.

Object definition and attributes#

A PointCloud is a Vector is a vector of 2D point geometries associated to numeric values from a main data column, and can also contain auxiliary data columns.

It inherits the main Vector attribute ds containing the geodataframe, and adds another main attribute data_column that identifies the name of the main data associated to the point geometries.

Additionally, new attributes such as point_count and new methods specific to point clouds are detailed further below.

Generic vector attributes and methods are inherited through the Vector object, such as bounds, crs, reproject() and crop().

Tip

The complete list of Vector attributes and methods can be found in the Vector section of the API.

Open and save#

A PointCloud is opened by instantiating the class with a str, a pathlib.Path, a geopandas.GeoDataFrame, a geopandas.GeoSeries or a shapely.Geometry.

import geoutils
import numpy as np

# Instantiate a point cloud from a filename on disk
filename_dem = geoutils.examples.get_path("coromandel_lidar")
pc = geoutils.PointCloud(filename_dem, data_column="Z")
pc

PointCloud(
  ds=              Z                         geometry
       0       816.078  POINT (1838933.295 5887999.962)
       1       818.823  POINT (1838930.054 5887999.966)
       2       811.059  POINT (1838933.316 5887999.967)
       3       818.151  POINT (1838928.638 5887999.961)
       4       817.599  POINT (1838926.098 5887999.995)
       ...         ...                              ...
       584343  827.438  POINT (1838913.663 5888000.872)
       584344  825.099  POINT (1838912.806 5888000.715)
       584345  825.419  POINT (1838912.889 5888001.284)
       584346  825.204   POINT (1838912.775 5888000.72)
       584347  822.421  POINT (1838911.847 5888000.565)
       
       [584348 rows x 2 columns]
  crs=COMPD_CS["NZGD2000 / New Zealand Transverse Mercator 2000 + NZVD2016 height",PROJCS["NZGD2000 / New Zealand Transverse Mercator 2000",GEOGCS["NZGD2000",DATUM["New_Zealand_Geodetic_Datum_2000",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6167"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4167"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",173],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",1600000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Northing",NORTH],AXIS["Easting",EAST],AUTHORITY["EPSG","2193"]],VERT_CS["NZVD2016 height",VERT_DATUM["New Zealand Vertical Datum 2016",2005,AUTHORITY["EPSG","1169"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Gravity-related height",UP],AUTHORITY["EPSG","7839"]]]
  bounds=BoundingBox(left=np.float64(1838792.524), bottom=np.float64(5887910.586), right=np.float64(1838937.061), top=np.float64(5888036.092)))

Create from arrays or tuples#

A PointCloud is created from three 1D arrays, from a Nx3 or 3xN array, or from an iterable of 3-tuples by calling the class methods from_xyz(), from_array() or from_tuples() respectively.

# From three 1D arrays
pc1 = geoutils.PointCloud.from_xyz(x=np.array([1, 4]), y=np.array([2, 5]), z=np.array([3, 6]), crs=4326)
# From a Nx3 array
pc2 = geoutils.PointCloud.from_array(np.array([[1, 2, 3], [4, 5, 6]]), crs=4326)
# From a iterable of 3-tuples
pc3 = geoutils.PointCloud.from_tuples([(1, 2, 3), (4, 5, 6)], crs=4326)

Gridding to a raster#

Gridding a PointCloud into a specific grid can be done through the grid() function, that applies a gridding scheme according to a given methods (inverse-distance weighting, delauney triangle interpolation).

# Grid the point cloud on a 100x100 grid on its extent
coords = (np.linspace(pc.bounds.left, pc.bounds.right, 100), np.linspace(pc.bounds.bottom, pc.bounds.top, 100))
rst = pc.grid(grid_coords=coords)

Arithmetic#

A PointCloud can be applied any pythonic arithmetic operation (+, -, /, //, *, **, %) with another PointCloud, ndarray or number. It will output one or two PointClouds. NumPy coercion rules apply for dtype. The operation is applied to the data_column of the point cloud.

# Add 1 and divide point cloud by 2
(pc1 + 1)/2

PointCloud(
  ds=     z     geometry
       0  2.0  POINT (1 2)
       1  3.5  POINT (4 5)
  crs=EPSG:4326
  bounds=BoundingBox(left=np.float64(1.0), bottom=np.float64(2.0), right=np.float64(4.0), top=np.float64(5.0)))

A PointCloud can also be applied any pythonic logical comparison operation (==, !=, >=, >, <=, <) with another PointCloud, ndarray or number. It will cast to a boolean PointCloud.

# What are point cloud pixels are larger than 20?
pc1 > 20

PointCloud(
  ds=       z     geometry
       0  False  POINT (1 2)
       1  False  POINT (4 5)
  crs=EPSG:4326
  bounds=BoundingBox(left=np.float64(1.0), bottom=np.float64(2.0), right=np.float64(4.0), top=np.float64(5.0)))

See Support of pythonic operators for more details.

Array interface#

A PointCloud can be applied any NumPy universal functions and most mathematical, logical or masked-array functions with another PointCloud, ndarray or number. The operation is applied to the data_column of the point cloud.

# Compute the element-wise square-root
np.sqrt(pc1)

PointCloud(
  ds=          z     geometry
       0  1.732051  POINT (1 2)
       1  2.449490  POINT (4 5)
  crs=EPSG:4326
  bounds=BoundingBox(left=np.float64(1.0), bottom=np.float64(2.0), right=np.float64(4.0), top=np.float64(5.0)))

Logical comparison functions will cast to a boolean PointCloud.

# Is the pointcloud close to another one within tolerance?

np.isclose(pc1, pc1+0.05, atol=0.1)

PointCloud(
  ds=      z     geometry
       0  True  POINT (1 2)
       1  True  POINT (4 5)
  crs=EPSG:4326
  bounds=BoundingBox(left=np.float64(1.0), bottom=np.float64(2.0), right=np.float64(4.0), top=np.float64(5.0)))

See Masked-array NumPy interface for more details.

Statistics#

Statistics of a point cloud, optionally subsetting to an inlier mask, can be computed using get_stats().

# Get mean, max and STD of the point cloud
pc.get_stats(["mean", "max", "std"])

{'mean': np.float64(819.6583450461027),
 'max': np.float64(849.053),
 'std': np.float64(16.355949189397542)}

A point cloud can also be quickly subsampled using subsample(), which considers only valid values, and returns either a point cloud or an array:

# Get 500 random points in the point cloud
pc_sub = pc.subsample(500)

See Statistics for more details.