6

Working with Python 3.6 (Anaconda) under Windows using a laspy fork that works with Python 3 (https://github.com/sethrh/laspy).

I am loading a LAS file to a GeoPandas GeoDataFrame to eventually perform a spatial join (using sjoin) with a polygon layer. I am sure there is a more elegant way to do what I am doing using the kludge below:

import numpy as np
from laspy.file import File
from pandas import DataFrame
from geopandas import GeoDataFrame
from shapely.geometry import Point

#Read LAS file
inFile = File("s428_7568.las", mode = "r")

#Import LAS into numpy array (X=raw integer value x=scaled float value)
lidar_points = np.array((inFile.x,inFile.y,inFile.z,inFile.intensity,
               inFile.raw_classification,inFile.scan_angle_rank)).transpose()

#Transform to pandas DataFrame
lidar_df=DataFrame(lidar_points)

#Transform to geopandas GeoDataFrame
crs = None
geometry = [Point(xyz) for xyz in zip(inFile.x,inFile.y,inFile.z)]
lidar_geodf = GeoDataFrame(lidar_df, crs=crs, geometry=geometry)
lidar_geodf.crs = {'init' :'epsg:2959'} # set correct spatial reference

So, is there a better way, more efficient, more pythonesque way of injecting the .las file points (from laspy) to the GeoPandas dataframe without passing through a numpy array?

0

2 Answers 2

5

So, is there a better way, more efficient, more pythonesque way of injecting the .las file points (from laspy) to the GeoPandas dataframe without passing through a numpy array?

No, and I'm not sure why you think this is the least efficient way to go. laspy is underneath the covers making a memoryview to the data and the Numpy array is a wrapper over that. There isn't much copying going on at all. You could provide more background on why you think this is the bottleneck.

Note that the latest 1.5.0 laspy release has included Seth's and other's patches to be compatible with Python 3.5/3.6.

2
  • I am a new to programming in Python so just wanted to ensure I was starting off correctly. Thanks for your insight! And yes, I saw on Git about version 1.5 which is great news as well. Apr 10, 2017 at 18:44
  • 1
    I am not sure that the list comprehension turning xyz coordinates to points is really efficient. This seems to be a bootleneck in my current, similar code.
    – CharlesG
    Oct 20, 2021 at 11:59
1

Not sure if this method was available in older versions of Geopandas, but this part

#Transform to geopandas GeoDataFrame
crs = None
geometry = [Point(xyz) for xyz in zip(inFile.x,inFile.y,inFile.z)]
lidar_geodf = GeoDataFrame(lidar_df, crs=crs, geometry=geometry)
lidar_geodf.crs = {'init' :'epsg:2959'} # set correct spatial reference

could be replaced with the more efficient:

from geopandas import GeoDataFrame, points_from_xy

#Transform to geopandas GeoDataFrame
lidar_geodf = GeoDataFrame(lidar_df,  geometry=points_from_xy(inFile.x, inFile.y, inFile.z, crs='epsg:2959'))

the function points_from_xy seems much more efficient than the list comprehension used

1
  • It was not available then, thank you for your post! Oct 19, 2023 at 14:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.