QGIS API Documentation  3.21.0-Master (909859188c)
qgseptpointcloudindex.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgspointcloudindex.cpp
3  --------------------
4  begin : October 2020
5  copyright : (C) 2020 by Peter Petrik
6  email : zilolv at gmail dot com
7  ***************************************************************************/
8 
9 /***************************************************************************
10  * *
11  * This program is free software; you can redistribute it and/or modify *
12  * it under the terms of the GNU General Public License as published by *
13  * the Free Software Foundation; either version 2 of the License, or *
14  * (at your option) any later version. *
15  * *
16  ***************************************************************************/
17 
18 #include "qgseptpointcloudindex.h"
19 #include <QFile>
20 #include <QFileInfo>
21 #include <QDir>
22 #include <QJsonArray>
23 #include <QJsonDocument>
24 #include <QJsonObject>
25 #include <QTime>
26 #include <QtDebug>
27 #include <QQueue>
28 
29 #include "qgseptdecoder.h"
31 #include "qgspointcloudrequest.h"
32 #include "qgspointcloudattribute.h"
33 #include "qgslogger.h"
34 #include "qgsfeedback.h"
35 #include "qgsmessagelog.h"
36 
38 
39 #define PROVIDER_KEY QStringLiteral( "ept" )
40 #define PROVIDER_DESCRIPTION QStringLiteral( "EPT point cloud provider" )
41 
42 QgsEptPointCloudIndex::QgsEptPointCloudIndex() = default;
43 
44 QgsEptPointCloudIndex::~QgsEptPointCloudIndex() = default;
45 
46 void QgsEptPointCloudIndex::load( const QString &fileName )
47 {
48  QFile f( fileName );
49  if ( !f.open( QIODevice::ReadOnly ) )
50  {
51  QgsMessageLog::logMessage( tr( "Unable to open %1 for reading" ).arg( fileName ) );
52  mIsValid = false;
53  return;
54  }
55 
56  const QDir directory = QFileInfo( fileName ).absoluteDir();
57  mDirectory = directory.absolutePath();
58 
59  const QByteArray dataJson = f.readAll();
60  bool success = loadSchema( dataJson );
61 
62  if ( success )
63  {
64  // try to import the metadata too!
65  QFile manifestFile( mDirectory + QStringLiteral( "/ept-sources/manifest.json" ) );
66  if ( manifestFile.open( QIODevice::ReadOnly ) )
67  {
68  const QByteArray manifestJson = manifestFile.readAll();
69  loadManifest( manifestJson );
70  }
71  }
72 
73  if ( success )
74  {
75  success = loadHierarchy();
76  }
77 
78  mIsValid = success;
79 }
80 
81 void QgsEptPointCloudIndex::loadManifest( const QByteArray &manifestJson )
82 {
83  QJsonParseError err;
84  // try to import the metadata too!
85  const QJsonDocument manifestDoc = QJsonDocument::fromJson( manifestJson, &err );
86  if ( err.error == QJsonParseError::NoError )
87  {
88  const QJsonArray manifestArray = manifestDoc.array();
89  // TODO how to handle multiple?
90  if ( ! manifestArray.empty() )
91  {
92  const QJsonObject sourceObject = manifestArray.at( 0 ).toObject();
93  const QString metadataPath = sourceObject.value( QStringLiteral( "metadataPath" ) ).toString();
94  QFile metadataFile( mDirectory + QStringLiteral( "/ept-sources/" ) + metadataPath );
95  if ( metadataFile.open( QIODevice::ReadOnly ) )
96  {
97  const QByteArray metadataJson = metadataFile.readAll();
98  const QJsonDocument metadataDoc = QJsonDocument::fromJson( metadataJson, &err );
99  if ( err.error == QJsonParseError::NoError )
100  {
101  const QJsonObject metadataObject = metadataDoc.object().value( QStringLiteral( "metadata" ) ).toObject();
102  if ( !metadataObject.empty() )
103  {
104  const QJsonObject sourceMetadata = metadataObject.constBegin().value().toObject();
105  mOriginalMetadata = sourceMetadata.toVariantMap();
106  }
107  }
108  }
109  }
110  }
111 }
112 
113 bool QgsEptPointCloudIndex::loadSchema( const QByteArray &dataJson )
114 {
115  QJsonParseError err;
116  const QJsonDocument doc = QJsonDocument::fromJson( dataJson, &err );
117  if ( err.error != QJsonParseError::NoError )
118  return false;
119  const QJsonObject result = doc.object();
120  mDataType = result.value( QLatin1String( "dataType" ) ).toString(); // "binary" or "laszip"
121  if ( mDataType != QLatin1String( "laszip" ) && mDataType != QLatin1String( "binary" ) && mDataType != QLatin1String( "zstandard" ) )
122  return false;
123 
124  const QString hierarchyType = result.value( QLatin1String( "hierarchyType" ) ).toString(); // "json" or "gzip"
125  if ( hierarchyType != QLatin1String( "json" ) )
126  return false;
127 
128  mSpan = result.value( QLatin1String( "span" ) ).toInt();
129  mPointCount = result.value( QLatin1String( "points" ) ).toDouble();
130 
131  // WKT
132  const QJsonObject srs = result.value( QLatin1String( "srs" ) ).toObject();
133  mWkt = srs.value( QLatin1String( "wkt" ) ).toString();
134 
135  // rectangular
136  const QJsonArray bounds = result.value( QLatin1String( "bounds" ) ).toArray();
137  if ( bounds.size() != 6 )
138  return false;
139 
140  const QJsonArray boundsConforming = result.value( QLatin1String( "boundsConforming" ) ).toArray();
141  if ( boundsConforming.size() != 6 )
142  return false;
143  mExtent.set( boundsConforming[0].toDouble(), boundsConforming[1].toDouble(),
144  boundsConforming[3].toDouble(), boundsConforming[4].toDouble() );
145  mZMin = boundsConforming[2].toDouble();
146  mZMax = boundsConforming[5].toDouble();
147 
148  const QJsonArray schemaArray = result.value( QLatin1String( "schema" ) ).toArray();
150 
151  for ( const QJsonValue &schemaItem : schemaArray )
152  {
153  const QJsonObject schemaObj = schemaItem.toObject();
154  const QString name = schemaObj.value( QLatin1String( "name" ) ).toString();
155  const QString type = schemaObj.value( QLatin1String( "type" ) ).toString();
156 
157  const int size = schemaObj.value( QLatin1String( "size" ) ).toInt();
158 
159  if ( type == QLatin1String( "float" ) && ( size == 4 ) )
160  {
162  }
163  else if ( type == QLatin1String( "float" ) && ( size == 8 ) )
164  {
166  }
167  else if ( size == 1 )
168  {
170  }
171  else if ( type == QLatin1String( "unsigned" ) && size == 2 )
172  {
174  }
175  else if ( size == 2 )
176  {
178  }
179  else if ( size == 4 )
180  {
182  }
183  else
184  {
185  // unknown attribute type
186  return false;
187  }
188 
189  double scale = 1.f;
190  if ( schemaObj.contains( QLatin1String( "scale" ) ) )
191  scale = schemaObj.value( QLatin1String( "scale" ) ).toDouble();
192 
193  double offset = 0.f;
194  if ( schemaObj.contains( QLatin1String( "offset" ) ) )
195  offset = schemaObj.value( QLatin1String( "offset" ) ).toDouble();
196 
197  if ( name == QLatin1String( "X" ) )
198  {
199  mOffset.set( offset, mOffset.y(), mOffset.z() );
200  mScale.set( scale, mScale.y(), mScale.z() );
201  }
202  else if ( name == QLatin1String( "Y" ) )
203  {
204  mOffset.set( mOffset.x(), offset, mOffset.z() );
205  mScale.set( mScale.x(), scale, mScale.z() );
206  }
207  else if ( name == QLatin1String( "Z" ) )
208  {
209  mOffset.set( mOffset.x(), mOffset.y(), offset );
210  mScale.set( mScale.x(), mScale.y(), scale );
211  }
212 
213  // store any metadata stats which are present for the attribute
214  AttributeStatistics stats;
215  if ( schemaObj.contains( QLatin1String( "count" ) ) )
216  stats.count = schemaObj.value( QLatin1String( "count" ) ).toInt();
217  if ( schemaObj.contains( QLatin1String( "minimum" ) ) )
218  stats.minimum = schemaObj.value( QLatin1String( "minimum" ) ).toDouble();
219  if ( schemaObj.contains( QLatin1String( "maximum" ) ) )
220  stats.maximum = schemaObj.value( QLatin1String( "maximum" ) ).toDouble();
221  if ( schemaObj.contains( QLatin1String( "count" ) ) )
222  stats.mean = schemaObj.value( QLatin1String( "mean" ) ).toDouble();
223  if ( schemaObj.contains( QLatin1String( "stddev" ) ) )
224  stats.stDev = schemaObj.value( QLatin1String( "stddev" ) ).toDouble();
225  if ( schemaObj.contains( QLatin1String( "variance" ) ) )
226  stats.variance = schemaObj.value( QLatin1String( "variance" ) ).toDouble();
227  mMetadataStats.insert( name, stats );
228 
229  if ( schemaObj.contains( QLatin1String( "counts" ) ) )
230  {
231  QMap< int, int > classCounts;
232  const QJsonArray counts = schemaObj.value( QLatin1String( "counts" ) ).toArray();
233  for ( const QJsonValue &count : counts )
234  {
235  const QJsonObject countObj = count.toObject();
236  classCounts.insert( countObj.value( QLatin1String( "value" ) ).toInt(), countObj.value( QLatin1String( "count" ) ).toInt() );
237  }
238  mAttributeClasses.insert( name, classCounts );
239  }
240  }
241  setAttributes( attributes );
242 
243  // save mRootBounds
244 
245  // bounds (cube - octree volume)
246  const double xmin = bounds[0].toDouble();
247  const double ymin = bounds[1].toDouble();
248  const double zmin = bounds[2].toDouble();
249  const double xmax = bounds[3].toDouble();
250  const double ymax = bounds[4].toDouble();
251  const double zmax = bounds[5].toDouble();
252 
253  mRootBounds = QgsPointCloudDataBounds(
254  ( xmin - mOffset.x() ) / mScale.x(),
255  ( ymin - mOffset.y() ) / mScale.y(),
256  ( zmin - mOffset.z() ) / mScale.z(),
257  ( xmax - mOffset.x() ) / mScale.x(),
258  ( ymax - mOffset.y() ) / mScale.y(),
259  ( zmax - mOffset.z() ) / mScale.z()
260  );
261 
262 
263 #ifdef QGIS_DEBUG
264  double dx = xmax - xmin, dy = ymax - ymin, dz = zmax - zmin;
265  QgsDebugMsgLevel( QStringLiteral( "lvl0 node size in CRS units: %1 %2 %3" ).arg( dx ).arg( dy ).arg( dz ), 2 ); // all dims should be the same
266  QgsDebugMsgLevel( QStringLiteral( "res at lvl0 %1" ).arg( dx / mSpan ), 2 );
267  QgsDebugMsgLevel( QStringLiteral( "res at lvl1 %1" ).arg( dx / mSpan / 2 ), 2 );
268  QgsDebugMsgLevel( QStringLiteral( "res at lvl2 %1 with node size %2" ).arg( dx / mSpan / 4 ).arg( dx / 4 ), 2 );
269 #endif
270 
271  return true;
272 }
273 
274 QgsPointCloudBlock *QgsEptPointCloudIndex::nodeData( const IndexedPointCloudNode &n, const QgsPointCloudRequest &request )
275 {
276  mHierarchyMutex.lock();
277  const bool found = mHierarchy.contains( n );
278  mHierarchyMutex.unlock();
279  if ( !found )
280  return nullptr;
281 
282  if ( mDataType == QLatin1String( "binary" ) )
283  {
284  const QString filename = QStringLiteral( "%1/ept-data/%2.bin" ).arg( mDirectory, n.toString() );
285  return QgsEptDecoder::decompressBinary( filename, attributes(), request.attributes(), scale(), offset() );
286  }
287  else if ( mDataType == QLatin1String( "zstandard" ) )
288  {
289  const QString filename = QStringLiteral( "%1/ept-data/%2.zst" ).arg( mDirectory, n.toString() );
290  return QgsEptDecoder::decompressZStandard( filename, attributes(), request.attributes(), scale(), offset() );
291  }
292  else if ( mDataType == QLatin1String( "laszip" ) )
293  {
294  const QString filename = QStringLiteral( "%1/ept-data/%2.laz" ).arg( mDirectory, n.toString() );
295  return QgsEptDecoder::decompressLaz( filename, attributes(), request.attributes(), scale(), offset() );
296  }
297  else
298  {
299  return nullptr; // unsupported
300  }
301 }
302 
303 QgsPointCloudBlockRequest *QgsEptPointCloudIndex::asyncNodeData( const IndexedPointCloudNode &n, const QgsPointCloudRequest &request )
304 {
305  Q_UNUSED( n );
306  Q_UNUSED( request );
307  Q_ASSERT( false );
308  return nullptr; // unsupported
309 }
310 
312 {
314 }
315 
316 qint64 QgsEptPointCloudIndex::pointCount() const
317 {
318  return mPointCount;
319 }
320 
321 QVariant QgsEptPointCloudIndex::metadataStatistic( const QString &attribute, QgsStatisticalSummary::Statistic statistic ) const
322 {
323  if ( !mMetadataStats.contains( attribute ) )
324  return QVariant();
325 
326  const AttributeStatistics &stats = mMetadataStats[ attribute ];
327  switch ( statistic )
328  {
330  return stats.count >= 0 ? QVariant( stats.count ) : QVariant();
331 
333  return std::isnan( stats.mean ) ? QVariant() : QVariant( stats.mean );
334 
336  return std::isnan( stats.stDev ) ? QVariant() : QVariant( stats.stDev );
337 
339  return stats.minimum;
340 
342  return stats.maximum;
343 
345  return stats.minimum.isValid() && stats.maximum.isValid() ? QVariant( stats.maximum.toDouble() - stats.minimum.toDouble() ) : QVariant();
346 
360  return QVariant();
361  }
362  return QVariant();
363 }
364 
365 QVariantList QgsEptPointCloudIndex::metadataClasses( const QString &attribute ) const
366 {
367  QVariantList classes;
368  const QMap< int, int > values = mAttributeClasses.value( attribute );
369  for ( auto it = values.constBegin(); it != values.constEnd(); ++it )
370  {
371  classes << it.key();
372  }
373  return classes;
374 }
375 
376 QVariant QgsEptPointCloudIndex::metadataClassStatistic( const QString &attribute, const QVariant &value, QgsStatisticalSummary::Statistic statistic ) const
377 {
378  if ( statistic != QgsStatisticalSummary::Count )
379  return QVariant();
380 
381  const QMap< int, int > values = mAttributeClasses.value( attribute );
382  if ( !values.contains( value.toInt() ) )
383  return QVariant();
384  return values.value( value.toInt() );
385 }
386 
387 bool QgsEptPointCloudIndex::loadHierarchy()
388 {
389  QQueue<QString> queue;
390  queue.enqueue( QStringLiteral( "0-0-0-0" ) );
391  while ( !queue.isEmpty() )
392  {
393  const QString filename = QStringLiteral( "%1/ept-hierarchy/%2.json" ).arg( mDirectory, queue.dequeue() );
394  QFile fH( filename );
395  if ( !fH.open( QIODevice::ReadOnly ) )
396  {
397  QgsDebugMsgLevel( QStringLiteral( "unable to read hierarchy from file %1" ).arg( filename ), 2 );
398  return false;
399  }
400 
401  const QByteArray dataJsonH = fH.readAll();
402  QJsonParseError errH;
403  const QJsonDocument docH = QJsonDocument::fromJson( dataJsonH, &errH );
404  if ( errH.error != QJsonParseError::NoError )
405  {
406  QgsDebugMsgLevel( QStringLiteral( "QJsonParseError when reading hierarchy from file %1" ).arg( filename ), 2 );
407  return false;
408  }
409 
410  const QJsonObject rootHObj = docH.object();
411  for ( auto it = rootHObj.constBegin(); it != rootHObj.constEnd(); ++it )
412  {
413  const QString nodeIdStr = it.key();
414  const int nodePointCount = it.value().toInt();
415  if ( nodePointCount < 0 )
416  {
417  queue.enqueue( nodeIdStr );
418  }
419  else
420  {
421  const IndexedPointCloudNode nodeId = IndexedPointCloudNode::fromString( nodeIdStr );
422  mHierarchyMutex.lock();
423  mHierarchy[nodeId] = nodePointCount;
424  mHierarchyMutex.unlock();
425  }
426  }
427  }
428  return true;
429 }
430 
431 bool QgsEptPointCloudIndex::isValid() const
432 {
433  return mIsValid;
434 }
435 
Represents a indexed point cloud node in octree.
static IndexedPointCloudNode fromString(const QString &str)
Creates node from string.
QString toString() const
Encode node to string.
This class represents a coordinate reference system (CRS).
static QgsCoordinateReferenceSystem fromWkt(const QString &wkt)
Creates a CRS from a WKT spatial ref sys definition string.
static void logMessage(const QString &message, const QString &tag=QString(), Qgis::MessageLevel level=Qgis::MessageLevel::Warning, bool notifyUser=true)
Adds a message to the log instance (and creates it if necessary).
Collection of point cloud attributes.
void push_back(const QgsPointCloudAttribute &attribute)
Adds extra attribute.
Attribute for point cloud data pair of name and size in bytes.
@ UShort
Unsigned short int 2 bytes.
Base class for handling loading QgsPointCloudBlock asynchronously.
Base class for storing raw data from point cloud nodes.
Represents packaged data bounds.
Point cloud data request.
QgsPointCloudAttributeCollection attributes() const
Returns attributes.
Statistic
Enumeration of flags that specify statistics to be calculated.
@ InterQuartileRange
Inter quartile range (IQR)
@ Median
Median of values.
@ Variety
Variety (count of distinct) values.
@ First
First value (since QGIS 3.6)
@ Last
Last value (since QGIS 3.6)
@ Minority
Minority of values.
@ ThirdQuartile
Third quartile.
@ Majority
Majority of values.
@ CountMissing
Number of missing (null) values.
@ Range
Range of values (max - min)
@ StDevSample
Sample standard deviation of values.
@ FirstQuartile
First quartile.
@ StDev
Standard deviation of values.
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:39
const QgsCoordinateReferenceSystem & crs