Merge branch 'get_rid_of_liblas' into 'master'

Get rid of liblas

See merge request Oslandia/py3dtiles!78

docs/install.rst

@@ -1,6 +1,20 @@
 Install
 -------
 
+Dependencies:
+- PDAL > 1.7
+
+From pypi
+~~~~~~~~~~~~
+
+`py3dtiles` is published on pypi.org.
+
+```
+pip install py3dtiles
+```
+
+
+
 From sources
 ~~~~~~~~~~~~
 

py3dtiles/convert.py

@@ -8,11 +8,12 @@
 from collections import namedtuple
 import pickle
 import zmq
-import pyproj
+from pyproj import CRS, Transformer
 import psutil
 import struct
 import concurrent.futures
 import argparse
+from py3dtiles.utils import SrsInMissingException
 from py3dtiles.points.transformations import rotation_matrix, angle_between_vectors, vector_product, inverse_matrix, scale_matrix, translation_matrix
 from py3dtiles.points.utils import compute_spacing, name_to_filename
 from py3dtiles.points.node import Node
@@ -27,11 +28,7 @@
 total_memory_MB = int(psutil.virtual_memory().total / (1024 * 1024))
 
 
-class SrsInMissingException(Exception):
-    pass
-
-
-def write_tileset(in_folder, out_folder, octree_metadata, offset, scale, projection, rotation_matrix, include_rgb):
+def write_tileset(in_folder, out_folder, octree_metadata, offset, scale, rotation_matrix, include_rgb):
     # compute tile transform matrix
     if rotation_matrix is None:
         transform = np.identity(4)
@@ -98,7 +95,7 @@ def make_rotation_matrix(z1, z2):
 OctreeMetadata = namedtuple('OctreeMetadata', ['aabb', 'spacing', 'scale'])
 
 
-def zmq_process(activity_graph, projection, node_store, octree_metadata, folder, write_rgb, verbosity):
+def zmq_process(activity_graph, transformer, node_store, octree_metadata, folder, write_rgb, verbosity):
     context = zmq.Context()
 
     # Socket to receive messages on
@@ -144,7 +141,7 @@ def zmq_process(activity_graph, projection, node_store, octree_metadata, folder,
                 command['offset_scale'],
                 command['portion'],
                 skt,
-                projection,
+                transformer,
                 verbosity)
         elif command[0] == b'pnts':
             command_type = 3
@@ -368,34 +365,30 @@ def convert(files,
     # read all input files headers and determine the aabb/spacing
     _, ext = os.path.splitext(files[0])
     init_reader_fn = las_reader.init if ext == '.las' else xyz_reader.init
-    infos = init_reader_fn(files, color_scale=color_scale, srs_in=srs_in)
+    infos = init_reader_fn(files, color_scale=color_scale, srs_in=srs_in, srs_out=srs_out)
 
     avg_min = infos['avg_min']
     rotation_matrix = None
     # srs stuff
-    projection = None
+    transformer = None
     if srs_out is not None:
-        p2 = pyproj.Proj(init='epsg:{}'.format(srs_out))
+        crs_out = CRS('epsg:{}'.format(srs_out))
         if srs_in is not None:
-            p1 = pyproj.Proj(init='epsg:{}'.format(srs_in))
-        else:
-            p1 = infos['srs_in']
-        if srs_in is None:
+            crs_in = CRS('epsg:{}'.format(srs_in))
+        elif infos['srs_in'] is None:
             raise SrsInMissingException('No SRS informations in the provided files')
-        projection = [p1, p2]
+        else:
+            crs_in = CRS(infos['srs_in'])
+        transformer = Transformer.from_crs(crs_in, crs_out)
 
-        bl = np.array(list(pyproj.transform(
-            projection[0], projection[1],
+        bl = np.array(list(transformer.transform(
             infos['aabb'][0][0], infos['aabb'][0][1], infos['aabb'][0][2])))
-        tr = np.array(list(pyproj.transform(
-            projection[0], projection[1],
+        tr = np.array(list(transformer.transform(
             infos['aabb'][1][0], infos['aabb'][1][1], infos['aabb'][1][2])))
-        br = np.array(list(pyproj.transform(
-            projection[0], projection[1],
+        br = np.array(list(transformer.transform(
             infos['aabb'][1][0], infos['aabb'][0][1], infos['aabb'][0][2])))
 
-        avg_min = np.array(list(pyproj.transform(
-            projection[0], projection[1],
+        avg_min = np.array(list(transformer.transform(
             avg_min[0], avg_min[1], avg_min[2])))
 
         x_axis = br - bl
@@ -501,7 +494,7 @@ def add_tasks_to_process(state, name, task, point_count):
     zmq_processes = [multiprocessing.Process(
         target=zmq_process,
         args=(
-            graph, projection, node_store, octree_metadata, outfolder, rgb, verbose)) for i in range(jobs)]
+            graph, transformer, node_store, octree_metadata, outfolder, rgb, verbose)) for i in range(jobs)]
 
     for p in zmq_processes:
         p.start()
@@ -658,7 +651,6 @@ def add_tasks_to_process(state, name, task, point_count):
                               octree_metadata,
                               avg_min,
                               root_scale,
-                              projection,
                               rotation_matrix,
                               rgb)
                 shutil.rmtree(working_dir)

py3dtiles/points/task/las_reader.py

@@ -1,13 +1,15 @@
+import json
 import numpy as np
 import math
 import traceback
 import laspy
-import pyproj
 import struct
+import subprocess
 from laspy.file import File
-import liblas
 from pickle import dumps as pdumps
 
+from py3dtiles.utils import SrsInMissingException
+
 
 def init(files, color_scale=None, srs_in=None, srs_out=None, fraction=100):
     aabb = None
@@ -51,12 +53,10 @@ def init(files, color_scale=None, srs_in=None, srs_out=None, fraction=100):
             pointcloud_file_portions += [(filename, p)]
 
         if (srs_out is not None and srs_in is None):
-            f = liblas.file.File(filename)
-            if (f.header.srs.proj4 is not None
-                    and f.header.srs.proj4 != ''):
-                srs_in = pyproj.Proj(f.header.srs.proj4)
-            else:
-                raise Exception('\'{}\' file doesn\'t contain srs information. Please use the --srs_in option to declare it.'.format(filename))
+            summary = json.loads(subprocess.check_output(['pdal', 'info', '--summary', filename]))['summary']
+            if 'srs' not in summary:
+                raise SrsInMissingException('\'{}\' file doesn\'t contain srs information. Please use the --srs_in option to declare it.'.format(filename))
+            srs_in = summary['srs']['proj4']
 
     return {
         'portions': pointcloud_file_portions,
@@ -68,7 +68,7 @@ def init(files, color_scale=None, srs_in=None, srs_out=None, fraction=100):
     }
 
 
-def run(_id, filename, offset_scale, portion, queue, projection, verbose):
+def run(_id, filename, offset_scale, portion, queue, transformer, verbose):
     '''
     Reads points from a las file
     '''
@@ -106,8 +106,8 @@ def run(_id, filename, offset_scale, portion, queue, projection, verbose):
             y = Y[start_offset:start_offset + num] * f.header.scale[1] + f.header.offset[1]
             z = Z[start_offset:start_offset + num] * f.header.scale[2] + f.header.offset[2]
 
-            if projection:
-                x, y, z = pyproj.transform(projection[0], projection[1], x, y, z)
+            if transformer:
+                x, y, z = transformer.transform(x, y, z)
 
             x = (x + offset_scale[0][0]) * offset_scale[1][0]
             y = (y + offset_scale[0][1]) * offset_scale[1][1]

py3dtiles/points/task/xyz_reader.py

@@ -1,7 +1,6 @@
 import numpy as np
 import math
 import traceback
-import pyproj
 import struct
 from pickle import dumps as pdumps
 
@@ -81,7 +80,7 @@ def init(files, color_scale=None, srs_in=None, srs_out=None, fraction=100):
     }
 
 
-def run(_id, filename, offset_scale, portion, queue, projection, verbose):
+def run(_id, filename, offset_scale, portion, queue, transformer, verbose):
     """
     Reads points from a xyz file
 
@@ -124,8 +123,8 @@ def run(_id, filename, offset_scale, portion, queue, projection, verbose):
 
             x, y, z = [points[:, c] for c in [0, 1, 2]]
 
-            if projection:
-                x, y, z = pyproj.transform(projection[0], projection[1], x, y, z)
+            if transformer:
+                x, y, z = transformer.transform(x, y, z)
 
             x = (x + offset_scale[0][0]) * offset_scale[1][0]
             y = (y + offset_scale[0][1]) * offset_scale[1][1]

py3dtiles/utils.py

@@ -1,15 +1,20 @@
 # -*- coding: utf-8 -*-
 
 import numpy as np
-import pyproj
+from pyproj import CRS, Transformer
 from .pnts import Pnts
 from .b3dm import B3dm
 
 
+class SrsInMissingException(Exception):
+    pass
+
+
 def convert_to_ecef(x, y, z, epsg_input):
-    inp = pyproj.Proj(init='epsg:{0}'.format(epsg_input))
-    outp = pyproj.Proj(init='epsg:4978')  # ECEF
-    return pyproj.transform(inp, outp, x, y, z)
+    inp = CRS('epsg:{0}'.format(epsg_input))
+    outp = CRS('epsg:4978')  # ECEF
+    transformer = Transformer.from_crs(inp, outp)
+    return transformer.transform(x, y, z)
 
 
 class TileContentReader(object):

setup.py

@@ -11,16 +11,15 @@
     'cython',
     'triangle',
     'psycopg2-binary',
-    'liblas',
     'laspy',
     'numba',
-    'pyproj',
     'psutil',
     'lz4',
     'pyzmq'
 )
 
 dev_requirements = (
+    'flake8',
     'pytest',
     'pytest-cov',
     'pytest-benchmark',

tests/test_convert.py

@@ -1,12 +1,21 @@
 # -*- coding: utf-8 -*-
 import os
-from pytest import approx, raises
+from pytest import approx, raises, fixture
 import shutil
 
 from py3dtiles import convert_to_ecef
 from py3dtiles.convert import convert, SrsInMissingException
 
 
+fixtures_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'fixtures')
+
+
+@fixture()
+def tmp_dir():
+    yield './tmp'
+    shutil.rmtree('./tmp', ignore_errors=True)
+
+
 def test_convert_to_ecef():
     # results tested with gdaltransform
     [x, y, z] = convert_to_ecef(-75.61200462622627,
@@ -18,25 +27,31 @@ def test_convert_to_ecef():
     approx(z, 4083122.83975827)
 
 
-def test_convert():
+def test_convert(tmp_dir):
     # just
-    convert(os.path.join(os.path.dirname(os.path.abspath(__file__)), './ripple.las'), outfolder='./tmp')
-    assert os.path.exists(os.path.join('tmp', 'tileset.json'))
-    assert os.path.exists(os.path.join('tmp', 'r0.pnts'))
-    shutil.rmtree('./tmp')
+    convert(os.path.join(os.path.dirname(os.path.abspath(__file__)), './ripple.las'), outfolder=tmp_dir)
+    assert os.path.exists(os.path.join(tmp_dir, 'tileset.json'))
+    assert os.path.exists(os.path.join(tmp_dir, 'r0.pnts'))
 
 
-def test_convert_without_srs():
+def test_convert_without_srs(tmp_dir):
     with raises(SrsInMissingException):
-        convert(os.path.join(os.path.dirname(os.path.abspath(__file__)), 'fixtures', 'without_srs.las'),
-                outfolder='./tmp',
+        convert(os.path.join(fixtures_dir, 'without_srs.las'),
+                outfolder=tmp_dir,
                 srs_out='4978')
-    assert not os.path.exists(os.path.join('tmp'))
+    assert not os.path.exists(os.path.join(tmp_dir))
 
-    convert(os.path.join(os.path.dirname(os.path.abspath(__file__)), 'fixtures', 'without_srs.las'),
-            outfolder='./tmp',
+    convert(os.path.join(fixtures_dir, 'without_srs.las'),
+            outfolder=tmp_dir,
             srs_in='3949',
             srs_out='4978')
-    assert os.path.exists(os.path.join('tmp', 'tileset.json'))
-    assert os.path.exists(os.path.join('tmp', 'r.pnts'))
-    shutil.rmtree('./tmp')
+    assert os.path.exists(os.path.join(tmp_dir, 'tileset.json'))
+    assert os.path.exists(os.path.join(tmp_dir, 'r.pnts'))
+
+
+def test_convert_with_srs(tmp_dir):
+    convert(os.path.join(fixtures_dir, 'with_srs.las'),
+            outfolder=tmp_dir,
+            srs_out='4978')
+    assert os.path.exists(os.path.join(tmp_dir, 'tileset.json'))
+    assert os.path.exists(os.path.join(tmp_dir, 'r.pnts'))