Source code for flusstools.lidartools.laspy_processor

from .laspy_config import *

[docs]class LasPoint: """Las file container to convert datasets to ESRI point shapefiles and/or GeoTIFFs. Args: las_file_name (str): Directory to and name of a las file. epsg (int): Authority Code - Geodetic Parameter Dataset ID (default: ``3857``). overwrite (bool): Overwrite existing shapefiles and/or GeoTIFFs (default: ``True``). use_attributes (str): Attributes (properties) to use from the las-file available in pattr ( (default: ``use_attributes="aciw"``). Attributes: las_file (laspy.file.File): A laspy file object attributes (str): Defined with ``use_attributes`` epsg (int): Authority code gdf (geopandas.GeoDataFrame): geopandas data frame containing all points of the las file with the properties (columns) defined by ``use_attributes`` offset (laspy.file.File().header.offset): Offset of las points (auto-read) overwrite (bool): Enable or disable overwriting existing files (default: ``True``) scale (laspy.file.File().header.scale): Scale of las points relative to the offset (auto-read) shapefile_name (str): The name and dicrectorty of a point shapefile where all las-file data is stored srs (osr.SpatialReference): The geo-spatial reference imported from ``epsg`` """ def __init__(self, las_file_name, epsg=3857, use_attributes="aciw", overwrite=True): self.las_file = laspy.file.File(las_file_name, mode="r") self.attributes = use_attributes self.epsg = epsg self.gdf = geopandas.GeoDataFrame() # void initialization self.offset = np.array(self.las_file.header.offset, dtype=np.float64) self.overwrite = overwrite self.scale = np.array(self.las_file.header.scale, dtype=np.float64) self.shapefile_name = "" self.srs = osr.SpatialReference() self.srs.ImportFromEPSG(epsg)"Using EPSG = %04i" % epsg) def __del__(self): self.las_file.close() del self.las_file def __repr__(self): return "%s" % self.__class__.__name__
[docs] def create_dem(self, target_file_name="", pixel_size=1.0, **kwargs): """Creates a digital elevation model (DEM) in GeoTIFF format from the *las* file points. Args: target_file_name (str): A file name including an existing directory where the dem will be created< must end on ``.tif``. pixel_size (float): The size of one pixel relative to the spatial reference system Keyword Args: src_shp_file_name (str): Name of a shapefile from which elevation information is to be extracted (default: name of the las-point shapefile) elevation_field_name (str): Name of the field from which elevation data is to be extracted (default: ``"elevation"``) interpolate_gap_pixels (bool): Fill empty pixels that are not touched by a shapefile point with interpolated values (default: ``False``) radius1 (float): Define the x-radius for interpolating pixels (default: ``-1``, corresponding to infinity). Only applicable ``with interpolate_gap_pixels``. radius2 (float): Define the y-radius for interpolating pixels (default: ``-1``, corresponding to infinity). Only applicable ``with interpolate_gap_pixels``. power (float): Power of the function for interpolating pixel values (default: ``1.0``, corresponding to linear). smoothing (float): Smoothing parameter for interpolating pixel values (default: ``0.0``). min_points (int): Minimum number of points to use for interpolation. If the interpolator cannot find at least ``min_points`` for a pixel, it assigns a ``no_data`` value to that pixel (default: ``0``). max_points (int): Maximum number of points to use for interpolation. The interpolator will not use more than ``max_points`` closest points to interpolate a pixel value (default: ``0``). Hint: This function works independently and does not require the prior creation of a shapefile. Returns: int: ``0`` if successful, otherwise ``-1`` """" * Creating GeoTIFF DEM %s ..." % target_file_name) default_keys = {"src_shp_file_name": self.shapefile_name, "elevation_field_name": "Elevation", "interpolate_gap_pixels": False, "radius1": -1, "radius2": -1, "power": 1.0, "smoothing": 0.0, "min_points": 0, "max_points": 0, } for k in default_keys.keys(): if kwargs.get(k): default_keys[k] = str(kwargs.get(k)) if not os.path.isfile(default_keys["src_shp_file_name"]):" * Need to create a point shapefile first (%s does not exist) ..." % default_keys["src_shp_file_name"]) self.export2shp(shapefile_name=default_keys["src_shp_file_name"]) if os.path.isfile(target_file_name) and self.overwrite:" -- Overwriting %s ..." % target_file_name) geo_utils.remove_tif(target_file_name) geo_utils.rasterize(default_keys["src_shp_file_name"], target_file_name, pixel_size=pixel_size, field_name=default_keys["elevation_field_name"], interpolate_gap_pixels=default_keys["interpolate_gap_pixels"], power=default_keys["power"], radius1=default_keys["radius1"], radius2=default_keys["radius2"], smoothing=default_keys["smoothing"], min_points=default_keys["min_points"], max_points=default_keys["max_points"])" -- Done.") return 0
[docs] def export2shp(self, **kwargs): """Converts las file points to a point shapefile. Keyword Args: shapefile_name (`str`): Optional shapefile name (must end on .shp). (default: ``'/this/dir/las_file_name.shp'``). Returns: str: ``/path/to/shapefile.shp``, which is a point shapefile created by the function. """ if kwargs.get("shapefile_name"): self.shapefile_name = kwargs.get("shapefile_name") else: self.shapefile_name = os.path.abspath("") + "/{0}.shp".format(self.las_file.filename) if os.path.isfile(self.shapefile_name) and self.overwrite is False:" * Using existing shapefile %s." % self.shapefile_name) return self.shapefile_name self._build_data_frame()" * Writing geopandas.GeoDataFrame to shapefile (%s) ..." % self.shapefile_name)" *** this action may take a while (0.5h per 1 million points)***") self.gdf.to_file(filename=self.shapefile_name, driver="ESRI Shapefile")" -- Done.") return self.shapefile_name
[docs] def get_file_info(self): """ Prints las file information to console.""" print("Point data formats in file:") for f in self.las_file.point_format: print(" -- %s" % print("File header info:") headers = [str( for spec in self.las_file.header.header_format] print(" -- " + ", ".join(headers))
def _build_data_frame(self): """ Builds the geopandas GeoDataFrame - auto-runs ``self._parse_attributes``.""" point_dict = self._parse_attributes() # for attr in self.pts_description: # if not"[x-z]", attr): # point_dict.update({attr: pts_df[attr]})" * Building geopandas.GeoDataFrame ...") self.gdf = geopandas.GeoDataFrame(pd.DataFrame(point_dict), crs="EPSG:%04i" % self.epsg)" -- Done.") def _get_xyz_array(self): """Extract x-y-z data from las records in a faster way than using ``las_file.x``, ``y``, or ``z``. Returns: ndarray: The DEM information extracted from the las file. """ pts = self.las_file.points['point'].copy().view(np.recarray) # read and transform data (from raw - fast than las_file.x) dem_array = np.empty((3, len(pts)), dtype=np.float64) dem_array[0, :] = pts.X * self.scale[0] + self.offset[0] dem_array[1, :] = pts.Y * self.scale[1] + self.offset[1] dem_array[2, :] = pts.Z * self.scale[2] + self.offset[2] return dem_array def _parse_attributes(self): """Parses attributes and append entries to point list."""" * Extracting transformed point coordinates ...") dem = self._get_xyz_array() point_dict = {"geometry": geopandas.points_from_xy(x=dem[0], y=dem[1], z=dem[2])} # add elevation field to facilitate DEM export point_dict.update({"elevation": self.las_file.z})" * Parsing and extracting user attributes of points ...") for attr in self.attributes: try: point_dict.update({wattr[attr]: self.las_file.__getattribute__(pattr[attr])})" -- added %s" % wattr[attr]) except AttributeError: logging.error("Non-existing attribute %s. Valid attributes are: %s" % (str(attr), dict2str(wattr))) except KeyError: logging.error("Non-existing las-file key %s - valid are: " % str(attr) + ", ".join(dir(self.las_file))) return point_dict