brycefrank/pyfor

Tools for analyzing aerial point clouds of forest data.

pyfor's Issues

Subversive pyproj error for Windows conda environment

Certain environment configurations in Windows have problematic pyproj installations that are unable to find local spatial reference files:

 File "_proj.pyx", line 129, in _proj.Proj.__cinit__
RuntimeError: b'no arguments in initialization list'

Oddly, this does not fail if ran from the command line directly, only from PyCharm

Cache miniconda travis directory

Due to the dependency stack, installation of the environment on travis is slow, it should be possible to cache arbitrary directories via .travis.yml for a speed up e.g. :

theochem/cgrid#8

Detected segments misprojected

Hello! I am getting a warning saying this: "This behavior has changed from < 0.3.1, points are now binned from the top left of the point cloud instead of the bottom right to cohere with arrays produced later". And my Shapefile result is not located in the place it should be while the projection set to it is correct. Do you know any solution to this?

Thanks for your help!

`create_index` should be multithreaded by default

pyfor/pyfor/collection.py

Lines 226 to 232 in 091a666

	def create_index(self):
	"""
	For each file in the collection, creates `.lax` files for spatial indexing using the default values.
	"""
	for las_path in self['las_path']:
	laxpy.file.init_lax(las_path)

Unnecessary axis inversion in DetectedTops plotting method

pyfor/pyfor/rasterizer.py

Line 291 in 67c24ac

fig.gca().invert_yaxis()

see #26

TIN option for surfaces

Hello again, apologies for the multiple posts.

Another option that would be very useful for us in our workflow is the possibility to use a TIN algorithm to create our surfaces, similar to this function in the lidR package (we are currently using).

https://www.rdocumentation.org/packages/lidR/versions/2.0.1/topics/p2r

Thank you!

Failing docs

RTD failed on latest master build due to excessive memory consumption.

Python 2

Is there any way to use pyfor with python2? I have to use it in a project including ArcPy which needs python2 to work.

Thank for the reply.

Point gridding can be simplified

The following gridding operation is too complex:

pyfor/pyfor/rasterizer.py

Lines 35 to 39 in 091a666

	bins_x = np.searchsorted(x_edges, self.cloud.data.points['x'], side='right') - 1
	bins_y = np.searchsorted(-y_edges, -self.cloud.data.points['y'], side='left', sorter=(-y_edges).argsort())
	self.cloud.data.points.loc[:, "bins_x"] = bins_x
	self.cloud.data.points.loc[:, "bins_y"] = bins_y

The above can be generalized to any arbitrary bin edges (irregular distances, etc), but is too complex for the case of regular bin edges. This can be replaced with a simple division by resolution and floor operation. The y axis needs to be flipped in some way first such that the bin indices are matrix style.

Summary output

It would be nice to have a summary output for Cloud and Grid objects for interactive data analysis. This would also force us to look at proper Cloud object creation. This would probably end up being a function/property of CloudData.

Something like:

cloud1 = cloud.Cloud("data/sample.las")
cloud1.las.summary

File Name: sample.las
File Size: 12 MB
Point Count: 1,123,456
Point Density: 15.2
etc

Writing a custom field -> python R portability.

i'm 99.9% this isn't a pyfor problem, but I saw some commits, so i thought i'd check with you first.

During my segmentation pipeline I make a new column "Tree" for individual tree ID. I like the visualization tools in the R lidR package, so i went go load it up there.

In python

pc.data.points.columns
Index(['x', 'y', 'z', 'intensity', 'return_num', 'classification', 'flag_byte',
       'scan_angle_rank', 'user_data', 'pt_src_id', 'bins_x', 'bins_y',
       'Tree'],
      dtype='object')

in R, we are missing that column. I'm trying to figure out if it got written.

> colnames(r@data)
 [1] "X"                          "Y"                          "Z"                         
 [4] "gpstime"                    "Intensity"                  "ReturnNumber"              
 [7] "NumberOfReturns"            "ScanDirectionFlag"          "EdgeOfFlightline"          
[10] "Classification"             "Synthetic_flag"             "Keypoint_flag"             
[13] "Withheld_flag"              "ScanAngle"                  "UserData"                  
[16] "PointSourceID"              "R"                          "G"                         
[19] "B"                          "reversible index (lastile)"

I'm pretty sure this is the lidR package dropping the column, but let me know if i'm missing something.

Your code for writing is pretty unambiguous.

pyfor/pyfor/cloud.py

Line 70 in 2108ac9

for dim in self.points:

So i think pyfor is safe there. I'm gonna checkin with lidR.

Add warning for large number of tiles for a retiling operation

Currently it is easy for a user to re-tile it too fine a resolution without warning. This leads to a bottleneck when computing area level geometries. A simple UserWarning in Retiler should suffice.

Travis fails on write .laz

======================================================================

ERROR: test_write (pyfortest.test_pyfor.LAZCloudTestCase)

----------------------------------------------------------------------

Traceback (most recent call last):

  File "/home/travis/build/brycefrank/pyfor/pyfortest/test_pyfor.py", line 171, in setUp

    self.test_cloud = cloud.Cloud(test_laz)

  File "/home/travis/build/brycefrank/pyfor/pyfor/cloud.py", line 88, in __init__

    las = laspy.file.File(self.filepath)

  File "/home/travis/miniconda/envs/pyfor_env/lib/python3.7/site-packages/laspy/file.py", line 64, in __init__

    self.open()

  File "/home/travis/miniconda/envs/pyfor_env/lib/python3.7/site-packages/laspy/file.py", line 75, in open

    self._reader = base.Reader(self.filename, mode=self._mode)

  File "/home/travis/miniconda/envs/pyfor_env/lib/python3.7/site-packages/laspy/base.py", line 273, in __init__

    self.setup_read_write(vlrs,evlrs, read_only=True)

  File "/home/travis/miniconda/envs/pyfor_env/lib/python3.7/site-packages/laspy/base.py", line 321, in setup_read_write

    self.data_provider.point_map()

  File "/home/travis/miniconda/envs/pyfor_env/lib/python3.7/site-packages/laspy/base.py", line 176, in point_map

    count=self.manager.header.point_records_count)

TypeError: a bytes-like object is required, not 'FakeMmap'

Interpolation error

Some distortion or error is present when using the interpolation function.

Here is a plotted array before interpolation using the code:

import matplotlib.pyplot as plt
import pyfor

A = pyfor.cloud.Cloud("sample.las")
A = a_grid.array("min", "z")
plt.imshow(A);
plt.colorbar()
plt.show()

Here is a plotted array after interpolation using the code:

A = a_grid.interpolate("z", "min")
plt.imshow(A);
plt.colorbar()
plt.show()

Area-based metrics. New deep learning tools.

Using this issue to do a bit of introductions. @brycefrank , meet @eayrey. Elias just published an awesome example of deep learning for area-based metrics. Bryce is working on a lidar forest python module. I've been pushing my friends at the Forest service to explore very similar ideas. Paper attached.
remotesensing-10-00649.pdf

If there is a trained model available (?) from

https://github.com/Eayrey/3D-Convolutional-Neural-Networks-with-LiDAR

it might be fairly straightforward to pass a PointCloud method to the trained model and get biomass/treecount predictions. Just an idea.

Warning for `force_extent` bbox

pyfor/pyfor/rasterizer.py

Line 188 in 1700e58

new_left, new_right, new_bot, new_top = bbox

It would be useful and perhaps cleaner to round the bbox to the nearest raster-compatible extent if
the user did not supply one. The user should be warned if this occurs.

Add heightbreak for all metrics

Right now metrics are computed using all points, with the exception of cover metrics. Would be useful to allow the user to define a global height break.

Eliminate `ResourceWarning`s from testing suite

e.g.:

ResourceWarning: Enable tracemalloc to get the object allocation traceback
/home/travis/miniconda/envs/pyfor_env/lib/python3.7/site-packages/llvmlite/ir/builder.py:910: ResourceWarning: unclosed file <_io.BufferedReader name='/home/travis/build/brycefrank/pyfor/pyfortest/data/mock_collection/1.las'>
  idx = [idx]

Probably a result of bad reading habits in Cloud

Creation of conda environment does not install rasterio

conda env create -f environment.yml

ResolvePackageNotFound:
  - rasterio[version='>=1.0.2']

`KrausPfeifer1998` classify fails

Traceback (most recent call last):
  File "C:\Users\bfrank70\Miniconda3\envs\pyfor_env\lib\site-packages\pandas\core\indexes\base.py", line 2897, in get_loc
    return self._engine.get_loc(key)
  File "pandas/_libs/index.pyx", line 107, in pandas._libs.index.IndexEngine.get_loc
  File "pandas/_libs/index.pyx", line 131, in pandas._libs.index.IndexEngine.get_loc
  File "pandas/_libs/hashtable_class_helper.pxi", line 1607, in pandas._libs.hashtable.PyObjectHashTable.get_item
  File "pandas/_libs/hashtable_class_helper.pxi", line 1614, in pandas._libs.hashtable.PyObjectHashTable.get_item
KeyError: 'v_i'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "C:\Users\bfrank70\Miniconda3\envs\pyfor_env\lib\unittest\case.py", line 59, in testPartExecutor
    yield
  File "C:\Users\bfrank70\Miniconda3\envs\pyfor_env\lib\unittest\case.py", line 605, in run
    testMethod()
  File "C:\Users\bfrank70\Programming\pyfor\pyfortest\test_pyfor.py", line 388, in test_classify
    self.test_kp_filter.classify(self.test_cloud)
  File "C:\Users\bfrank70\Programming\pyfor\pyfor\ground_filter.py", line 269, in classify
    grid.cloud.data.points["classification"][grid.cloud.data.points['v_i'] <= self.g + self.w] = ground_int
  File "C:\Users\bfrank70\Miniconda3\envs\pyfor_env\lib\site-packages\pandas\core\frame.py", line 2980, in __getitem__
    indexer = self.columns.get_loc(key)
  File "C:\Users\bfrank70\Miniconda3\envs\pyfor_env\lib\site-packages\pandas\core\indexes\base.py", line 2899, in get_loc
    return self._engine.get_loc(self._maybe_cast_indexer(key))
  File "pandas/_libs/index.pyx", line 107, in pandas._libs.index.IndexEngine.get_loc
  File "pandas/_libs/index.pyx", line 131, in pandas._libs.index.IndexEngine.get_loc
  File "pandas/_libs/hashtable_class_helper.pxi", line 1607, in pandas._libs.hashtable.PyObjectHashTable.get_item
  File "pandas/_libs/hashtable_class_helper.pxi", line 1614, in pandas._libs.hashtable.PyObjectHashTable.get_item
KeyError: 'v_i'

Scanning files for a collection should be generalized to `.las/z` files

pyfor/pyfor/collection.py

Lines 36 to 37 in 091a666

	las_path_init = [[os.path.join(root, file) for file in files] for root, dirs, files in os.walk(las_dir)][0]
	las_path_init = [las_path for las_path in las_path_init if las_path.endswith('.las')]

Error in median pit filtering for canopy height model.

Sorry to issue blitz.

same tile as #15

import pyfor
from shapely import geometry
from matplotlib import pyplot

def createPolygon(xmin,xmax,ymin,ymax):
    '''
    Convert a pandas row into a polygon bounding box
    ''' 
    
    p1 = geometry.Point(xmin,ymax)
    p2 = geometry.Point(xmax,ymax)
    p3 = geometry.Point(xmax,ymin)
    p4 = geometry.Point(xmin,ymin)
    
    pointList = [p1, p2, p3, p4, p1]
    
    poly = geometry.Polygon([[p.x, p.y] for p in pointList])
    
    return poly

# Grab a polygon
window_xmin=407690.2
window_xmax=407715.2
window_ymin=3291477.2
window_ymax=3291452.2

#Create shapely polygon
poly=createPolygon(window_xmin, window_xmax, window_ymax,window_ymin)

#Read lidar tile
pc=pyfor.cloud.Cloud("/Users/ben/Documents/DeepForest/data/NEON_D03_OSBS_DP1_407000_3291000_classified_point_cloud.laz")

#Filter
pc.filter(min = -5, max = 100, dim = "z")

#Clip to geographic extent
clipped=pc.clip(poly)

#Compute Canopy model
chm = clipped.chm(cell_size = 0.5, interp_method = "nearest")
chm.plot()

chm = clipped.chm(cell_size = 0.5, interp_method = "nearest", pit_filter = "median", kernel_size = 3)

chm.plot()

I could be wrong, but I don't think passing a 3X3 filter across the image could create this image. Some kind of artifact being created. Should I try pulling from a different branch?

Raster Class

Create a raster class that will serve as a wrapper for the hairier raster I/O functions of GDAL without having to use rasterio as a dependency. This class will represent a given manifestation of a grid object for some function (i.e. the min Z value in each cell) and will be the last step before output. This also will provide methods for importing raster data.

i.e.:

las1 = pyfor.cloud.Cloud("my.las")
las1.grid(1).metrics("max", "z").write_raster("mytif.tif") #.metrics in this case returns a Raster

or in the case of input

rast1 = pyfor.raster.Raster("my.tif")
rast1.plot() # plot a 2d raster in matplotlib
rast1.array # extract a numpy array of the raster

Merge update_docs.sh into Makefile

pyfor/docs/update_docs.sh

Lines 1 to 8 in 0b9c8a9

	#!/bin/bash
	# This script autoregenerates pyfor's documentation for updates to master.
	DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" >/dev/null && pwd )"
	sphinx-apidoc -f -o api/ $DIR/../pyfor --separate
	cd $DIR
	make html
	cd -

Point cloud colors

Hello pyfor builders! Thank you for the tool.

Any chance you can include RGB values when reading LAS/LAZ files?
or is there a specific reason why you eliminate them?
It should be a matter of adding 'red','green','blue' within the list named dims within the method _get_las_points.

Thank you for the consideration!

Watershed segmentation

Hi Bryce,

I am running your script of watershed segmentation but I discover that before save the shapefile, the shapefile is rotated. I defined the projection but it haven't improved it. Could you help me?

Thanks!
Carmen

Unclear argument parameter name in `retile_raster`

pyfor/pyfor/collection.py

Line 141 in 908f863

def retile_raster(self, cell_size, original_tile_size, buffer=0):

The argument original_tile_size is not completely accurate, it is more a of a "target" size for the retiling operation

Migrate to `conda-forge` testing

Currently, the pyfor conda-forge feedstock does not run the testing suite, as there were some issues with file locations. These errors should be resolved and testing for the master branch replaced with testing of the current conda-forge build. This implies changing the badge on the master branch,

Single Photon Lidar data

Hello,

First of all as a fellow PhD student working in forestry, thank you so much for your efforts on this package. I work on tree species identification at the individual tree level and have developped a bunch of python scripts to calculate 3D and intensity features from ALS data. Your package is exactly what we need to make our process as Python-centric as possible.

I am working with Leica SPL100 data (in LAS 1.4 format which I know causes problems in some open source implementations, since you are using laspy I believe this is not the problem), the concept of a scan angle does not exist for this type of sensor. When I try to load a LAS file with this data, I get:

LaspyException: Dimension: scan_angle_rank not found.

Instead of an exception would it be possible to just ignore it?

Once again many thanks for your work on this, much appreciated!

Flipping the canopy height model?

Hey Bryce,

A tiny question today, i'll really try not to lean too much on ya.

I'm working on deep learning on joint RGB and lidar imagery, so I need to crop and overlay both data at the same window. I'm all good to go, but it looks like the canopy height model is mirrored over the y axis when plotting.

When overlaying the RGB (image) and LIDAR (CHM), its easy to see.

        fig, ax = pyplot.subplots()
        ax.imshow(image[:,:,::-1])
        ax.matshow(CHM.array,alpha=0.4)
        pyplot.show()

Then flip the CHM over the y axis

        fig, ax = pyplot.subplots()
        ax.imshow(np.flip(image[:,:,::-1],0))
        ax.matshow(CHM.array,alpha=0.4)
        pyplot.show()

I think you already noticed this in pyfor:

pyfor/pyfor/rasterizer.py

Line 170 in eb158be

fig.gca().invert_yaxis()

Do you know what's happening here? Is it

1. That matplotlib "shows" the y axis inverted. So you are fixing it visually (this might be dangerous?).
2. The matrix itself is flipped, because of the way pyfor calculates the CHM, so you invert it for visualization
3. Neither 1 or 2, meaning my image is flipped (seems less likely).

It matters because the next step is to concatenate the CHM to the 3 channel rgb to create a 4 channel input to the deep net. I don't want to stick it on upside down.

normalization produces blank cloud (missing data?)

I've noticed that on occasion the normalization function can strip the entire cloud of points.

For example here it works.

import pyfor
from matplotlib import pyplot as plt
pc=pyfor.cloud.Cloud("/Users/ben/Documents/DeepLidar/data/SJER/SJER_006.laz")

pc.plot()
pc.normalize(3)
pc.plot()

But here it doesn't

import pyfor
from matplotlib import pyplot as plt
pc=pyfor.cloud.Cloud("/Users/ben/Documents/DeepLidar/data/SJER/SJER_004.laz")

pc.plot()
pc.normalize(3)
pc.plot()

These were crops from a larger tile (I can reprocess if I must). It looks like some missing data in the corner might be to blame? Just a guess. But then again

pc.data.z.isna().sum()
0

pc.data.points.z
0      NaN
1      NaN
2      NaN
3      NaN
4      NaN
5      NaN
6      NaN

data can be found here: https://github.com/weecology/NeonTreeEvaluation/tree/master/SJER

Documentation build failing

Documentation is not building for some reason:

https://readthedocs.org/projects/pyfor-pdal-u/builds/7077493/

Reported to this thread:

readthedocs/readthedocs.org#3996

Seems to be an issue with RTD itself, the last build of the documentation is still hosted.

Easier way to check units, useful for threshold in watershed

Hey bryce - this is looking really great! Amazing job. I just came back to it for a way to pass watershed clusters to my deep learning tree classifier (still in progress). One thing i noticed is that its not immediately easy to know the map units of the vertical dimension. It might be useful to have a summary method for the cloud object that prints some essentials.

Similar to

https://rdrr.io/cran/lidR/man/summary.LAS.html

That way I know what units to input to the threshold tolerance for watershed.

Can't write to file when acting as class object?

I ran into a strange error that might be demonstrative of something larger. Forgive the lack of reproducibility.

Consider this pyfor object. It was just clipped from a larger plot.

generator.clipped_las
<pyfor.cloud.Cloud object at 0x14e07c438>

It looks normal

generator.clipped_las.data.points.head()
     index           x            y   ...    bins_x  bins_y  bins_z
0  2336576  315006.189  4091000.072   ...         6      37       0
1  2336577  315005.268  4091000.146   ...         5      37       0
2  2336578  315004.345  4091000.222   ...         4      37       0
3  2336579  315003.436  4091000.298   ...         3      37       0
4  2336580  315002.527  4091000.370   ...         2      37       0

and I can see the method

generator.clipped_las.write
<bound method Cloud.write of <pyfor.cloud.Cloud object at 0x14e07c438>>

but then

generator.clipped_las.write(path="/Users/Ben/Desktop/test.laz")
Traceback (most recent call last):
  Debug Probe, prompt 99, line 1
  File "/Users/ben/miniconda3/envs/DeepLidar_dask/lib/python3.6/site-packages/pyfor/cloud.py", line 350, in write
    self.data.write(path)
builtins.AttributeError: 'CloudData' object has no attribute 'write'

Nothing is wrong with my writer

pc=pyfor.cloud.Cloud("/Users/ben/Documents/DeepLidar/data/TEAK/TEAK_044.laz")
pc.write("/Users/Ben/Desktop/test2.laz")

Nor is it the clipping. Here is the same behavior with the original tile.

generator.lidar_tile.write(path="/Users/Ben/Desktop/fulltest.laz")
Traceback (most recent call last):
  Debug Probe, prompt 107, line 1
  File "/Users/ben/miniconda3/envs/DeepLidar_dask/lib/python3.6/site-packages/pyfor/cloud.py", line 350, in write
    self.data.write(path)
builtins.AttributeError: 'CloudData' object has no attribute 'write'

some use of super class methods in pyfor? Something is weird here. Happy to dig.

Cryptic error when clipping a particular polygon and point set

bool llvm::APInt::operator==(const llvm::APInt&) const: Assertion `BitWidth == RHS.BitWidth && "Comparison requires equal bit widths"' failed.

Profile Canopy Height model

pyfor/pyfor/cloud.py

Line 304 in d2bb376

def chm(self, cell_size, interp_method=None, pit_filter=None, kernel_size=3):

I've been trying to increase performance of the canopy height model. I had assumed that it was the cell size resolution that led to slow performance on big data sets. This was incorrect.

Here is a toy example.
Here is my loading function, but it shouldn't matter.

def load_lidar(laz_path):
    try:
        pc=pyfor.cloud.Cloud(laz_path)
        pc.extension=".las"    
    except FileNotFoundError:
        print("Failed loading path: %s" %(laz_path))
        
    #normalize and filter
    zhang_filter=pyfor.ground_filter.Zhang2003(pc, cell_size=1)
    zhang_filter.normalize()    
    
    #Quick filter for unreasonable points.
    pc.filter(min = -5, max = pc.data.points.z.quantile(0.995), dim = "z")    
    
    #Check dim
    assert (not pc.data.points.shape[0] == 0), "Lidar tile is empty!"

    return pc

import cProfile, pstats
cp = cProfile.Profile()

from DeepForest import Lidar

cp.enable()
lidar_tile=Lidar.load_lidar("/Users/ben/Documents/DeepLidar/data/SJER/SJER_002.laz")
chm = lidar_tile.chm(cell_size = 0.1 , interp_method = "nearest", pit_filter = "median", kernel_size = 11)
cp.disable()

         5226 function calls (5150 primitive calls) in 0.064 seconds

   Ordered by: cumulative time, call count
   List reduced from 517 to 40 due to restriction <40>

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.000    0.000    0.064    0.064 cloud.py:304(chm)
        1    0.000    0.000    0.049    0.049 cloud.py:157(grid)
        1    0.000    0.000    0.049    0.049 rasterizer.py:18(__init__)
        2    0.000    0.000    0.047    0.024 frame.py:3105(__setitem__)
        2    0.000    0.000    0.047    0.023 frame.py:3182(_set_item)
        2    0.000    0.000    0.046    0.023 generic.py:2633(_check_setitem_copy)
        1    0.046    0.046    0.046    0.046 {built-in method gc.collect}
        1    0.000    0.000    0.014    0.014 rasterizer.py:48(raster)

Add the pit filter and interpolation.

import cProfile, pstats
cp = cProfile.Profile()

from DeepForest import Lidar

lidar_tile=Lidar.load_lidar("/Users/ben/Documents/DeepLidar/data/SJER/SJER_002.laz")

cp.enable()
chm = lidar_tile.chm(cell_size = 0.1 , interp_method = "nearest", pit_filter = "median", kernel_size = 11)
cp.disable()

This behavior has changed from < 0.3.1, points are now binned from the top left of the point cloud instead of the bottom right to cohere with arrays produced later.

         203631 function calls (201530 primitive calls) in 1.750 seconds

   Ordered by: cumulative time, call count
   List reduced from 1184 to 40 due to restriction <40>

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.000    0.000    1.753    1.753 cloud.py:304(chm)
        1    0.000    0.000    1.529    1.529 rasterizer.py:244(pit_filter)
    162/1    0.002    0.000    0.796    0.796 <frozen importlib._bootstrap>:966(_find_and_load)
    162/1    0.001    0.000    0.796    0.796 <frozen importlib._bootstrap>:936(_find_and_load_unlocked)
    153/1    0.002    0.000    0.796    0.796 <frozen importlib._bootstrap>:651(_load_unlocked)
    123/1    0.001    0.000    0.795    0.795 <frozen importlib._bootstrap_external>:672(exec_module)
    219/1    0.000    0.000    0.795    0.795 <frozen importlib._bootstrap>:211(_call_with_frames_removed)
    336/1    0.034    0.000    0.795    0.795 {built-in method builtins.exec}
        1    0.046    0.046    0.795    0.795 __init__.py:308(<module>)
        1    0.000    0.000    0.733    0.733 signaltools.py:854(medfilt)
        1    0.733    0.733    0.733    0.733 {built-in method scipy.signal.sigtools._order_filterND}
  890/344    0.001    0.000    0.661    0.002 <frozen importlib._bootstrap>:997(_handle_fromlist)
        1    0.000    0.000    0.351    0.351 filter_design.py:2(<module>)
        1    0.032    0.032    0.350    0.350 __init__.py:266(<module>)
        1    0.004    0.004    0.272    0.272 _peak_finding.py:3(<module>)
        1    0.020    0.020    0.267    0.267 __init__.py:342(<module>)
        1    0.030    0.030    0.243    0.243 _minimize.py:8(<module>)
        1    0.004    0.004    0.217    0.217 stats.py:156(<module>)

0.064 / 1.17 = 18x drop.

My real data is 360MB 1km tile. I can report general times for that, but i stopped. I think its non-linear, the cell size = 0.1 takes about 150 seconds without pit filter, so that would be 45 minutes. I waited that long and it was still going. If you give me a push I can fork and play around. My gut feeling is this isn't right, clearly there is a cost of passing a kernel over the array, but 18x feels off. What do you think is the bottleneck here?

Silent failure of cloud class? Laz support?

Hey Bryce, does pyfor still support .laz?

I just git cloned, and made a new conda env

MacBook-Pro:pyfor ben$ conda activate pyfor_env
(pyfor_env) MacBook-Pro:pyfor ben$ python
Python 3.6.6 | packaged by conda-forge | (default, Jul 26 2018, 09:55:02)
[GCC 4.2.1 Compatible Apple LLVM 6.1.0 (clang-602.0.53)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import pyfor
>>> pc=pyfor.cloud.Cloud("/Users/ben/Documents/DeepLidar/data/SJER/SJER_002.laz")
>>> pc

pc.data
Traceback (most recent call last):
File "", line 1, in
AttributeError: 'Cloud' object has no attribute 'data'

https://www.dropbox.com/s/i1dnfuvl5y4x06m/SJER_003.laz?dl=0

I've tried with a few files.

It all looks good in R. I think its well formatted.

> a<-readLAS("/Users/ben/Dropbox/Weecology/NEON/SJER_003.laz")
> a
class        : LAS (LASF v1.3)
memory       : 889.8 Kb 
extent       : 257388 , 257428 , 4111280 , 4111320 (xmin, xmax, ymin, ymax)
area         : 1596.612 m² (convex hull)
points       : 11241 points,  NA pulses
density      : 7.04 points/m²,  NA pulses/m²
field names  : X Y Z gpstime Intensity ReturnNumber NumberOfReturns ScanDirectionFlag EdgeOfFlightline Classification Synthetic_flag Keypoint_flag Withheld_flag ScanAngle UserData PointSourceID 
coord. ref.  : NA 

Tree Segmentation

Hi Bryce!

Im somewhat new to python so please pardon my lack of knowledge.
Im trying to do a crown segmentation of a park at my university as part of my master thesis.
When I perform the watershed segmentation with:

segmentation = chm.watershed_seg(min_distance=4, threshold_abs=4).plot()

I get:

/peak.py:185: RuntimeWarning: invalid value encountered in greater
  mask &= image > max(thresholds)

Also I dont get any polygons on the output. Its just the CHM again. See figure..

When I run:

segmentation.segments[segmentation.segments["raster_val"] > 0]

I get:

'NoneType' object has no attribute 'segments'

To be honest, Im a bit clueless whats actually wrong here.
If you could help me it would be really appreciated.

Thanks so much in advance!
David

Interpolation extent issue

Hello,

I am trying to understand the behaviour of the interpolation when I am creating surfaces such as a DTM. As mentionned in my email discussion, our point clouds are already classified so I filter out the ground returns then use the following to create the DTM raster.

chm = tile.chm(cell_size=0.25, interp_method = "linear")

We are processing 1600 buffered 1km^2 las tiles (40m buffer) to create DTM, DSM and CHM which then gets used by our ITC segmentation script. The rasters are at 25cm resolution so the output raster tiles should be 4160*4160. Our previous method using the lidR package in R produced this result.

When I use pyfor (and another python package called whitebox), the rasters have uneven dimensions such as 4079*4159 for example. I understand what is happening, the ground las points do not reach the corner of the declared extent of the las file and the interpolation stops a bit short. Here are some screenshots that explain the issue

LAS ground points (1000m * 1000m)

lidR output (4160*4160)

pyfor output (4079*4079)

whitebox output for comparison (4080*4080)

Your package does a better job as it uses the points at the top to interpolate whereas wbt does not

Arcmap TIN display

I believe what may be happening is the scipy interpolation is not getting the extent from the LAS file but calculating it from the actual points which makes it fall short of the 'official' extent. Do not know what R does different in this respect but this would be the desired output.

I will be digging into the code eventually to see if I can force this myself.

Thank you!

New version of Hdf5 causing crashes in travis

Warning! ***HDF5 library version mismatched error***
The HDF5 header files used to compile this application do not match
the version used by the HDF5 library to which this application is linked.
Data corruption or segmentation faults may occur if the application continues.
This can happen when an application was compiled by one version of HDF5 but
linked with a different version of static or shared HDF5 library.
You should recompile the application or check your shared library related
settings such as 'LD_LIBRARY_PATH'.
You can, at your own risk, disable this warning by setting the environment
variable 'HDF5_DISABLE_VERSION_CHECK' to a value of '1'.
Setting it to 2 or higher will suppress the warning messages totally.
Headers are 1.10.2, library is 1.10.1

Normalizing after clipping.

A bit hard to make reproducible, but give me a sense for whether this is avoidable behavior.

I can normalize a point cloud before, but not after clipping. My thought was that it would be much faster to normalize the clipped cloud, since its a dramatically smaller area. The process will need to be 10,000s of times, so performance is somewhat of an issue, but I can work around it if needed.

laz_path
'data/SJER/NEON_D17_SJER_DP1_257000_4110000_classified_point_cloud_colorized.laz'

Original data, non-normalized

pc.data.points.z.min()
58.948
pc.data.points.z.max()
630.018

Normalizes and crops just fine

pc.normalize(5)
pc.filter(min = -5, max = 100, dim = "z")
clipped=pc.clip(poly)
clipped.data.points.z.min()
-0.7969999999999686
clipped.data.points.z.max()
11.634999999999991

But you can't first crop and then normalize

clipped=pc.clip(poly)

clipped.normalize(1)

Traceback (most recent call last):
  Debug Probe, prompt 131, line 3
    '''
  File "/Users/ben/miniconda3/envs/retinanet/lib/python3.6/site-packages/pyfor/cloud.py", line 262, in normalize
    filter.normalize(cell_size)
  File "/Users/ben/miniconda3/envs/retinanet/lib/python3.6/site-packages/pyfor/ground_filter.py", line 243, in normalize
    bem = self.bem(cell_size)
  File "/Users/ben/miniconda3/envs/retinanet/lib/python3.6/site-packages/pyfor/ground_filter.py", line 217, in bem
    ground_cloud = self.ground_points
  File "/Users/ben/miniconda3/envs/retinanet/lib/python3.6/site-packages/pyfor/ground_filter.py", line 206, in ground_points
    ground = self._filter()
  File "/Users/ben/miniconda3/envs/retinanet/lib/python3.6/site-packages/pyfor/ground_filter.py", line 196, in _filter
    return self.cloud.data.points.loc[ix[ground_bins]]
  File "/Users/ben/miniconda3/envs/retinanet/lib/python3.6/site-packages/pandas/core/indexing.py", line 1478, in __getitem__
    return self._getitem_axis(maybe_callable, axis=axis)
  File "/Users/ben/miniconda3/envs/retinanet/lib/python3.6/site-packages/pandas/core/indexing.py", line 1901, in _getitem_axis
    return self._getitem_iterable(key, axis=axis)
  File "/Users/ben/miniconda3/envs/retinanet/lib/python3.6/site-packages/pandas/core/indexing.py", line 1143, in _getitem_iterable
    self._validate_read_indexer(key, indexer, axis)
  File "/Users/ben/miniconda3/envs/retinanet/lib/python3.6/site-packages/pandas/core/indexing.py", line 1206, in _validate_read_indexer
    key=key, axis=self.obj._get_axis_name(axis)))
builtins.KeyError: 'None of [[35506496. 41295890. 35505383. ... 44116653. 44116654. 44116655.]] are in the [index]'

Maybe need to reset index for pandas frame?

IndexError when creating a VoxelGrid

IndexError                                Traceback (most recent call last)
<ipython-input-10-d572ce203027> in <module>()
     10         chm = cld.chm(cell_size = 1.64, interp_method = "nearest", pit_filter = "median")
---> 11         voxel_raster = voxel_grid.voxel_raster("count", "z")
     12         convex_hull = cld.convex_hull

~/miniconda3/envs/pyfor_env/lib/python3.6/site-packages/pyfor/voxelizer.py in voxel_raster(self, func, dim)
     43         # Set the values of the grid
---> 44         voxel_grid[cells["bins_x"], cells["bins_y"], cells["bins_z"]] = cells[dim]
     45 

IndexError: index 49 is out of bounds for axis 1 with size 49

Filetype check does not work with capitalized file extension

pyfor will error on, for example

pc = pyfor.cloud.Cloud("test.LAS")

Implement single source for package version

e.g. https://packaging.python.org/guides/single-sourcing-package-version/

Update axis labels on plot of cropped cloud.

I think the plot function needs to re-calculate the extent of the cloud.las.points objects before plotting. If you crop and then plot, the same extent is used?

Tile is here:
https://github.com/weecology/TreeSegmentation/blob/master/data/training/NEON_D03_OSBS_DP1_407000_3291000_classified_point_cloud.laz

import pyfor
from shapely import geometry

def createPolygon(xmin,xmax,ymin,ymax):
    '''
    Convert a pandas row into a polygon bounding box
    ''' 
    
    p1 = geometry.Point(xmin,ymax)
    p2 = geometry.Point(xmax,ymax)
    p3 = geometry.Point(xmax,ymin)
    p4 = geometry.Point(xmin,ymin)
    
    pointList = [p1, p2, p3, p4, p1]
    
    poly = geometry.Polygon([[p.x, p.y] for p in pointList])
    
    return poly

#Grab some coords.
window_xmin=407690.2
window_xmax=407715.2
window_ymin=3291477.2
window_ymax=3291452.2

#Create shapely polygon
poly=createPolygon(window_xmin, window_xmax, window_ymin, window_ymax)

#Read lidar tile
pc=pyfor.cloud.Cloud("/Users/ben/Documents/DeepForest/data/NEON_D03_OSBS_DP1_407000_3291000_classified_point_cloud.laz")

pc.filter(min = -5, max = 100, dim = "z")

#View result

#Clip to geographic extent
clipped=pc.clip(poly)

clipped.plot()

Version

(pyfor_env) MacBook-Pro:pyfor ben$ git log
commit a893ddb23d9192ea4f4b49d32ef1583ed6446dbe (HEAD -> master, origin/master, origin/HEAD)
Author: brycefrank <[email protected]>
Date:   Sat Jun 23 11:16:15 2018 -0700

    test data fix

Clip multiple features from single Cloud

Hi Bryce,

I would like to clip my point cloud with a shapefile with several polygons. It is possible to do with polyclip? How can I save this point cloud clipped?

Thanks!
Carmen

Filtering points

I'm going to move our convo here for others to be able see. There will be a few of these as I get up and running.

Use case. A tile has some outlier points I want to remove (z < 0) and (z < 100). What's the best way to filter based on cloud.las arrays. in lidR this would be analogous to

lasfilter

tile<-tile %>% lasfilter(Z < 100)
tile<-tile %>% lasfilter(Z > 0)

I can implement this if you give me a bit of a push.

We would need to subset the numpy array based on the condition. Then recalculate the cloud object to have the header update. Anything else?

Travis YML branch specification has no effect

pyfor/.travis.yml

Lines 3 to 6 in 7ac8604

	only:
	- master
	- 0.2.1
	- 0.3.0

Branches other than these build on travis.

Add avenues for creating bare earth models for already-classified points

This could easily be added directly to Zhang2003 and KrausPfeifer1998, should be accompanied with additions to the user manual.

Pandas error on plot after clipping

After clipping a point cloud and plotting that clipped point cloud the following warning appears:

A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy