QGIS API Documentation  3.18.1-Zürich (202f1bf7e5)
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  bool success = loadSchema( f );
59  if ( success )
60  {
61  success = loadHierarchy();
62  }
63 
64  mIsValid = success;
65 }
66 
67 bool QgsEptPointCloudIndex::loadSchema( QFile &f )
68 {
69  const QByteArray dataJson = f.readAll();
70  QJsonParseError err;
71  const QJsonDocument doc = QJsonDocument::fromJson( dataJson, &err );
72  if ( err.error != QJsonParseError::NoError )
73  return false;
74  const QJsonObject result = doc.object();
75  mDataType = result.value( QLatin1String( "dataType" ) ).toString(); // "binary" or "laszip"
76  if ( mDataType != QLatin1String( "laszip" ) && mDataType != QLatin1String( "binary" ) && mDataType != QLatin1String( "zstandard" ) )
77  return false;
78 
79  const QString hierarchyType = result.value( QLatin1String( "hierarchyType" ) ).toString(); // "json" or "gzip"
80  if ( hierarchyType != QLatin1String( "json" ) )
81  return false;
82 
83  mSpan = result.value( QLatin1String( "span" ) ).toInt();
84  mPointCount = result.value( QLatin1String( "points" ) ).toInt();
85 
86  // WKT
87  const QJsonObject srs = result.value( QLatin1String( "srs" ) ).toObject();
88  mWkt = srs.value( QLatin1String( "wkt" ) ).toString();
89 
90  // rectangular
91  const QJsonArray bounds = result.value( QLatin1String( "bounds" ) ).toArray();
92  if ( bounds.size() != 6 )
93  return false;
94 
95  const QJsonArray boundsConforming = result.value( QLatin1String( "boundsConforming" ) ).toArray();
96  if ( boundsConforming.size() != 6 )
97  return false;
98  mExtent.set( boundsConforming[0].toDouble(), boundsConforming[1].toDouble(),
99  boundsConforming[3].toDouble(), boundsConforming[4].toDouble() );
100  mZMin = boundsConforming[2].toDouble();
101  mZMax = boundsConforming[5].toDouble();
102 
103  const QJsonArray schemaArray = result.value( QLatin1String( "schema" ) ).toArray();
105 
106  for ( const QJsonValue &schemaItem : schemaArray )
107  {
108  const QJsonObject schemaObj = schemaItem.toObject();
109  const QString name = schemaObj.value( QLatin1String( "name" ) ).toString();
110  const QString type = schemaObj.value( QLatin1String( "type" ) ).toString();
111 
112  int size = schemaObj.value( QLatin1String( "size" ) ).toInt();
113 
114  if ( type == QLatin1String( "float" ) && ( size == 4 ) )
115  {
117  }
118  else if ( type == QLatin1String( "float" ) && ( size == 8 ) )
119  {
121  }
122  else if ( size == 1 )
123  {
125  }
126  else if ( type == QLatin1String( "unsigned" ) && size == 2 )
127  {
129  }
130  else if ( size == 2 )
131  {
133  }
134  else if ( size == 4 )
135  {
137  }
138  else
139  {
140  // unknown attribute type
141  return false;
142  }
143 
144  double scale = 1.f;
145  if ( schemaObj.contains( QLatin1String( "scale" ) ) )
146  scale = schemaObj.value( QLatin1String( "scale" ) ).toDouble();
147 
148  double offset = 0.f;
149  if ( schemaObj.contains( QLatin1String( "offset" ) ) )
150  offset = schemaObj.value( QLatin1String( "offset" ) ).toDouble();
151 
152  if ( name == QLatin1String( "X" ) )
153  {
154  mOffset.set( offset, mOffset.y(), mOffset.z() );
155  mScale.set( scale, mScale.y(), mScale.z() );
156  }
157  else if ( name == QLatin1String( "Y" ) )
158  {
159  mOffset.set( mOffset.x(), offset, mOffset.z() );
160  mScale.set( mScale.x(), scale, mScale.z() );
161  }
162  else if ( name == QLatin1String( "Z" ) )
163  {
164  mOffset.set( mOffset.x(), mOffset.y(), offset );
165  mScale.set( mScale.x(), mScale.y(), scale );
166  }
167 
168  // store any metadata stats which are present for the attribute
169  AttributeStatistics stats;
170  if ( schemaObj.contains( QLatin1String( "count" ) ) )
171  stats.count = schemaObj.value( QLatin1String( "count" ) ).toInt();
172  if ( schemaObj.contains( QLatin1String( "minimum" ) ) )
173  stats.minimum = schemaObj.value( QLatin1String( "minimum" ) ).toDouble();
174  if ( schemaObj.contains( QLatin1String( "maximum" ) ) )
175  stats.maximum = schemaObj.value( QLatin1String( "maximum" ) ).toDouble();
176  if ( schemaObj.contains( QLatin1String( "count" ) ) )
177  stats.mean = schemaObj.value( QLatin1String( "mean" ) ).toDouble();
178  if ( schemaObj.contains( QLatin1String( "stddev" ) ) )
179  stats.stDev = schemaObj.value( QLatin1String( "stddev" ) ).toDouble();
180  if ( schemaObj.contains( QLatin1String( "variance" ) ) )
181  stats.variance = schemaObj.value( QLatin1String( "variance" ) ).toDouble();
182  mMetadataStats.insert( name, stats );
183 
184  if ( schemaObj.contains( QLatin1String( "counts" ) ) )
185  {
186  QMap< int, int > classCounts;
187  const QJsonArray counts = schemaObj.value( QLatin1String( "counts" ) ).toArray();
188  for ( const QJsonValue &count : counts )
189  {
190  const QJsonObject countObj = count.toObject();
191  classCounts.insert( countObj.value( QLatin1String( "value" ) ).toInt(), countObj.value( QLatin1String( "count" ) ).toInt() );
192  }
193  mAttributeClasses.insert( name, classCounts );
194  }
195  }
196  setAttributes( attributes );
197 
198  // try to import the metadata too!
199 
200  QFile manifestFile( mDirectory + QStringLiteral( "/ept-sources/manifest.json" ) );
201  if ( manifestFile.open( QIODevice::ReadOnly ) )
202  {
203  const QByteArray manifestJson = manifestFile.readAll();
204  const QJsonDocument manifestDoc = QJsonDocument::fromJson( manifestJson, &err );
205  if ( err.error == QJsonParseError::NoError )
206  {
207  const QJsonArray manifestArray = manifestDoc.array();
208  // TODO how to handle multiple?
209  if ( ! manifestArray.empty() )
210  {
211  const QJsonObject sourceObject = manifestArray.at( 0 ).toObject();
212  const QString metadataPath = sourceObject.value( QStringLiteral( "metadataPath" ) ).toString();
213  QFile metadataFile( mDirectory + QStringLiteral( "/ept-sources/" ) + metadataPath );
214  if ( metadataFile.open( QIODevice::ReadOnly ) )
215  {
216  const QByteArray metadataJson = metadataFile.readAll();
217  const QJsonDocument metadataDoc = QJsonDocument::fromJson( metadataJson, &err );
218  if ( err.error == QJsonParseError::NoError )
219  {
220  const QJsonObject metadataObject = metadataDoc.object().value( QStringLiteral( "metadata" ) ).toObject();
221  if ( !metadataObject.empty() )
222  {
223  const QJsonObject sourceMetadata = metadataObject.constBegin().value().toObject();
224  mOriginalMetadata = sourceMetadata.toVariantMap();
225  }
226  }
227  }
228  }
229  }
230  }
231 
232  // save mRootBounds
233 
234  // bounds (cube - octree volume)
235  double xmin = bounds[0].toDouble();
236  double ymin = bounds[1].toDouble();
237  double zmin = bounds[2].toDouble();
238  double xmax = bounds[3].toDouble();
239  double ymax = bounds[4].toDouble();
240  double zmax = bounds[5].toDouble();
241 
242  mRootBounds = QgsPointCloudDataBounds(
243  ( xmin - mOffset.x() ) / mScale.x(),
244  ( ymin - mOffset.y() ) / mScale.y(),
245  ( zmin - mOffset.z() ) / mScale.z(),
246  ( xmax - mOffset.x() ) / mScale.x(),
247  ( ymax - mOffset.y() ) / mScale.y(),
248  ( zmax - mOffset.z() ) / mScale.z()
249  );
250 
251 
252 #ifdef QGIS_DEBUG
253  double dx = xmax - xmin, dy = ymax - ymin, dz = zmax - zmin;
254  QgsDebugMsgLevel( QStringLiteral( "lvl0 node size in CRS units: %1 %2 %3" ).arg( dx ).arg( dy ).arg( dz ), 2 ); // all dims should be the same
255  QgsDebugMsgLevel( QStringLiteral( "res at lvl0 %1" ).arg( dx / mSpan ), 2 );
256  QgsDebugMsgLevel( QStringLiteral( "res at lvl1 %1" ).arg( dx / mSpan / 2 ), 2 );
257  QgsDebugMsgLevel( QStringLiteral( "res at lvl2 %1 with node size %2" ).arg( dx / mSpan / 4 ).arg( dx / 4 ), 2 );
258 #endif
259 
260  return true;
261 }
262 
263 QgsPointCloudBlock *QgsEptPointCloudIndex::nodeData( const IndexedPointCloudNode &n, const QgsPointCloudRequest &request )
264 {
265  if ( !mHierarchy.contains( n ) )
266  return nullptr;
267 
268  if ( mDataType == QLatin1String( "binary" ) )
269  {
270  QString filename = QStringLiteral( "%1/ept-data/%2.bin" ).arg( mDirectory, n.toString() );
271  return QgsEptDecoder::decompressBinary( filename, attributes(), request.attributes() );
272  }
273  else if ( mDataType == QLatin1String( "zstandard" ) )
274  {
275  QString filename = QStringLiteral( "%1/ept-data/%2.zst" ).arg( mDirectory, n.toString() );
276  return QgsEptDecoder::decompressZStandard( filename, attributes(), request.attributes() );
277  }
278  else if ( mDataType == QLatin1String( "laszip" ) )
279  {
280  QString filename = QStringLiteral( "%1/ept-data/%2.laz" ).arg( mDirectory, n.toString() );
281  return QgsEptDecoder::decompressLaz( filename, attributes(), request.attributes() );
282  }
283  else
284  {
285  return nullptr; // unsupported
286  }
287 }
288 
290 {
292 }
293 
294 int QgsEptPointCloudIndex::pointCount() const
295 {
296  return mPointCount;
297 }
298 
299 QVariant QgsEptPointCloudIndex::metadataStatistic( const QString &attribute, QgsStatisticalSummary::Statistic statistic ) const
300 {
301  if ( !mMetadataStats.contains( attribute ) )
302  return QVariant();
303 
304  const AttributeStatistics &stats = mMetadataStats[ attribute ];
305  switch ( statistic )
306  {
308  return stats.count >= 0 ? QVariant( stats.count ) : QVariant();
309 
311  return std::isnan( stats.mean ) ? QVariant() : QVariant( stats.mean );
312 
314  return std::isnan( stats.stDev ) ? QVariant() : QVariant( stats.stDev );
315 
317  return stats.minimum;
318 
320  return stats.maximum;
321 
323  return stats.minimum.isValid() && stats.maximum.isValid() ? QVariant( stats.maximum.toDouble() - stats.minimum.toDouble() ) : QVariant();
324 
338  return QVariant();
339  }
340  return QVariant();
341 }
342 
343 QVariantList QgsEptPointCloudIndex::metadataClasses( const QString &attribute ) const
344 {
345  QVariantList classes;
346  const QMap< int, int > values = mAttributeClasses.value( attribute );
347  for ( auto it = values.constBegin(); it != values.constEnd(); ++it )
348  {
349  classes << it.key();
350  }
351  return classes;
352 }
353 
354 QVariant QgsEptPointCloudIndex::metadataClassStatistic( const QString &attribute, const QVariant &value, QgsStatisticalSummary::Statistic statistic ) const
355 {
356  if ( statistic != QgsStatisticalSummary::Count )
357  return QVariant();
358 
359  const QMap< int, int > values = mAttributeClasses.value( attribute );
360  if ( !values.contains( value.toInt() ) )
361  return QVariant();
362  return values.value( value.toInt() );
363 }
364 
365 bool QgsEptPointCloudIndex::loadHierarchy()
366 {
367  QQueue<QString> queue;
368  queue.enqueue( QStringLiteral( "0-0-0-0" ) );
369  while ( !queue.isEmpty() )
370  {
371  const QString filename = QStringLiteral( "%1/ept-hierarchy/%2.json" ).arg( mDirectory, queue.dequeue() );
372  QFile fH( filename );
373  if ( !fH.open( QIODevice::ReadOnly ) )
374  {
375  QgsDebugMsgLevel( QStringLiteral( "unable to read hierarchy from file %1" ).arg( filename ), 2 );
376  return false;
377  }
378 
379  QByteArray dataJsonH = fH.readAll();
380  QJsonParseError errH;
381  const QJsonDocument docH = QJsonDocument::fromJson( dataJsonH, &errH );
382  if ( errH.error != QJsonParseError::NoError )
383  {
384  QgsDebugMsgLevel( QStringLiteral( "QJsonParseError when reading hierarchy from file %1" ).arg( filename ), 2 );
385  return false;
386  }
387 
388  const QJsonObject rootHObj = docH.object();
389  for ( auto it = rootHObj.constBegin(); it != rootHObj.constEnd(); ++it )
390  {
391  QString nodeIdStr = it.key();
392  int nodePointCount = it.value().toInt();
393  if ( nodePointCount < 0 )
394  {
395  queue.enqueue( nodeIdStr );
396  }
397  else
398  {
400  mHierarchy[nodeId] = nodePointCount;
401  }
402  }
403  }
404  return true;
405 }
406 
407 bool QgsEptPointCloudIndex::isValid() const
408 {
409  return mIsValid;
410 }
411 
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::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 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