The georeferenced point cloud

The georeferenced point cloud#

In GeoUtils, georeferenced point clouds are exposed through the following interfaces:

For new workflows, we recommend using the accessors, and using either only accessors or only GeoUtils objects.

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

Important

The pc accessor on a Dataset supports processing several point clouds at once. See DataArray or Dataset accessors?.

Except for the differences outlined in Accessors pc versus PointCloud, the accessors and the PointCloud object provide the same geospatial operations, described in the feature pages linked in Feature and scalability overview and listed in Point cloud API.

Below, we describe how to open, create and save point clouds with the accessors, followed by the attributes and behaviour of the PointCloud object.

Accessors pc versus PointCloud#

The main differences are the following:

  1. A PointCloud applies arithmetic to its main data column, while the accessors require selecting the column or variable explicitly.

  2. A PointCloud supports Multiprocessing, while the accessors support Dask, allowing lazy graph building.

For numerical analysis, we recommend the Xarray accessor, and for geometric operations, the GeoPandas accessor. See Xarray or GeoPandas for points? for details.

Why a separate point object?#

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 supports 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 like LiDAR.

Point cloud accessors#

Opening files#

A point cloud is opened with open_pointcloud(), which returns a DataArray by default. Its values contain the main measurement, while auxiliary attributes are stored as coordinates along the same point dimension.

import geoutils

filename_dem = geoutils.examples.get_path("coromandel_lidar")

# Open elevation and classification as lazy numeric arrays from the LAS file
points = geoutils.open_pointcloud(filename_dem, data_name="Z", columns=["classification"], chunks=5000)
points

<xarray.DataArray 'Z' (point: 584348)> Size: 5MB
dask.array<concatenate, shape=(584348,), dtype=float64, chunksize=(5000,), chunktype=numpy.ndarray>
Coordinates:
    x               (point) float64 5MB dask.array<chunksize=(5000,), meta=np.ndarray>
    y               (point) float64 5MB dask.array<chunksize=(5000,), meta=np.ndarray>
    classification  (point) uint8 584kB dask.array<chunksize=(5000,), meta=np.ndarray>
Dimensions without coordinates: point
Attributes:
    crs:                COMPOUNDCRS["NZGD2000 / New Zealand Transverse Mercat...
    geometry_z:         False
    point_record_size:  30
    laz_chunk_size:     50000

Passing as_object_type="dataset" stores several measurements as separate variables, while as_object_type="geodataframe" returns a GeoDataFrame for working with point geometries. For LAS and GeoParquet files, columns selects the main attributes ("main"), all attributes ("all") or a list of names. The CRS, X/Y coordinates and main measurement are always read together.

By default, the Xarray data is read only when needed (load=False), without requiring Dask. We can select points with isel(point=...) before loading their values. This reads only the selected LAS/LAZ records, intersecting COPC nodes or intersecting GeoParquet row groups. Passing chunks= enables Dask execution and sets the number of points per chunk, while load=True reads the data into memory when opening the file.

For a GeoDataFrame, load=False uses Dask-GeoPandas, and chunks= sets the number of points per partition. Passing parse_sensor_metadata=True adds sensor information parsed from the filename.

Constructing from arrays#

A point cloud is created from X/Y coordinate and value arrays with from_xyz(), which also accepts auxiliary attributes:

points = geoutils.pointcloud.from_xyz([0, 1], [2, 3], [4, 5], crs=32631, data_name="height")

As with open_pointcloud(), passing as_object_type="dataset" or "geodataframe" returns the corresponding object type. Dask arrays remain lazy during construction.

Writing points#

A point cloud is saved to file with to_file(), to_las() or to_parquet(). For GeoParquet, a lazy array can be written in row groups of a chosen size or as a directory of partition files:

# Write lazy point arrays in groups of 5000 rows
points.pc.to_parquet("elevations.parquet", chunks=5000)

# Reopen the saved active measurement and auxiliary attributes
reopened = geoutils.open_pointcloud("elevations.parquet", columns="all", chunks=5000)

Writing runs the tasks needed to save the data, including for lazy arrays. The GeoParquet file stores the point coordinates, CRS and main measurement. See Formats and cloud storage for choosing formats and partition sizes, and Data objects and accessors for converting between array and geometry interfaces.

The PointCloud object#

Object definition and attributes#

A PointCloud is a vector of 2D point geometries associated with numeric values from a named data attribute, and can also contain auxiliary attributes.

It inherits the main Vector attribute gdf containing the geodataframe, and adds another main attribute data_name 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 bbox, 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_name="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"]]]
  bbox=BoundingBox(left=np.float64(1838792.524), bottom=np.float64(5887910.586), right=np.float64(1838937.061), top=np.float64(5888036.092)))

See Downsampling and overviews for details on reducing point cloud density when opening or plotting, including how GeoUtils selects deterministic point samples.

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)

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_name 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
  bbox=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
  bbox=BoundingBox(left=np.float64(1.0), bottom=np.float64(2.0), right=np.float64(4.0), top=np.float64(5.0)))

See Python 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_name 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
  bbox=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
  bbox=BoundingBox(left=np.float64(1.0), bottom=np.float64(2.0), right=np.float64(4.0), top=np.float64(5.0)))

See NumPy functions and reductions for more details.