diff --git a/python/core/auto_generated/pointcloud/qgspointcloudstatistics.sip.in b/python/core/auto_generated/pointcloud/qgspointcloudstatistics.sip.in
index 3464b3fbe96..5bcd99bb3b1 100644
--- a/python/core/auto_generated/pointcloud/qgspointcloudstatistics.sip.in
+++ b/python/core/auto_generated/pointcloud/qgspointcloudstatistics.sip.in
@@ -92,6 +92,16 @@ If no matching statistic is available then NaN will be returned.
     void combineWith( const QgsPointCloudStatistics &stats );
 %Docstring
 Merges the current statistics with the statistics from ``stats``
+%End
+
+    QByteArray toStatisticsJson() const;
+%Docstring
+Converts the current statistics object into JSON object
+%End
+
+    static QgsPointCloudStatistics fromStatisticsJson( QByteArray stats );
+%Docstring
+Creates a statistics object from the JSON object ``stats``
 %End
 
 };
diff --git a/src/core/pointcloud/qgscopcpointcloudindex.cpp b/src/core/pointcloud/qgscopcpointcloudindex.cpp
index 0357694f15a..d6870a3a48e 100644
--- a/src/core/pointcloud/qgscopcpointcloudindex.cpp
+++ b/src/core/pointcloud/qgscopcpointcloudindex.cpp
@@ -16,10 +16,14 @@
  ***************************************************************************/
 
 #include "qgscopcpointcloudindex.h"
+
+#include <fstream>
 #include <QFile>
 #include <QtDebug>
 #include <QQueue>
 #include <QMutexLocker>
+#include <QJsonDocument>
+#include <QJsonObject>
 
 #include "qgseptdecoder.h"
 #include "qgslazdecoder.h"
@@ -33,6 +37,7 @@
 
 #include "lazperf/lazperf.hpp"
 #include "lazperf/readers.hpp"
+#include "lazperf/vlr.hpp"
 
 ///@cond PRIVATE
 
@@ -181,6 +186,69 @@ bool QgsCopcPointCloudIndex::loadHierarchy()
   return true;
 }
 
+bool QgsCopcPointCloudIndex::writeStatistics( QgsPointCloudStatistics &stats )
+{
+  if ( mLazInfo->version() != qMakePair<uint8_t, uint8_t>( 1, 4 ) )
+  {
+    // EVLR isn't supported in the first place
+    QgsMessageLog::logMessage( tr( "Can't write statistics to \"%1\": laz version != 1.4" ).arg( mFileName ) );
+    return false;
+  }
+
+  QByteArray statisticsEvlrData = fetchCopcStatisticsEvlrData();
+  if ( !statisticsEvlrData.isEmpty() )
+  {
+    QgsMessageLog::logMessage( tr( "Can't write statistics to \"%1\": file already contains COPC statistics!" ).arg( mFileName ) );
+    return false;
+  }
+
+  lazperf::evlr_header statsEvlrHeader;
+  statsEvlrHeader.user_id = "qgis";
+  statsEvlrHeader.record_id = 0;
+  statsEvlrHeader.description = "Contains calculated statistics";
+  QByteArray statsJson = stats.toStatisticsJson();
+  statsEvlrHeader.data_length = statsJson.size();
+
+  // Save the EVLRs to the end of the original file (while erasing the exisitng EVLRs in the file)
+  mCopcFile.close();
+  std::fstream copcFile;
+  copcFile.open( QgsLazDecoder::toNativePath( mFileName ), std::ios_base::binary | std::iostream::in | std::iostream::out );
+  if ( copcFile.is_open() && copcFile.good() )
+  {
+    // Write the new number of EVLRs
+    lazperf::header14 header = mLazInfo->header();
+    header.evlr_count = header.evlr_count + 1;
+    copcFile.seekp( 0 );
+    header.write( copcFile );
+
+    // Append EVLR data to the end
+    copcFile.seekg( 0, std::ios::end );
+
+    statsEvlrHeader.write( copcFile );
+    copcFile.write( statsJson.data(), statsEvlrHeader.data_length );
+  }
+  else
+  {
+    QgsMessageLog::logMessage( tr( "Couldn't open COPC file \"%1\" to write statistics" ).arg( mFileName ) );
+    return false;
+  }
+  copcFile.close();
+  mCopcFile.open( QgsLazDecoder::toNativePath( mFileName ), std::ios::binary );
+  return true;
+}
+
+QgsPointCloudStatistics QgsCopcPointCloudIndex::readStatistics()
+{
+  QByteArray statisticsEvlrData = fetchCopcStatisticsEvlrData();
+
+  if ( statisticsEvlrData.isEmpty() )
+  {
+    return QgsPointCloudStatistics();
+  }
+
+  return QgsPointCloudStatistics::fromStatisticsJson( statisticsEvlrData );
+}
+
 bool QgsCopcPointCloudIndex::isValid() const
 {
   return mIsValid;
@@ -295,4 +363,33 @@ void QgsCopcPointCloudIndex::copyCommonProperties( QgsCopcPointCloudIndex *desti
   destination->mLazInfo.reset( new QgsLazInfo( *mLazInfo ) );
 }
 
+QByteArray QgsCopcPointCloudIndex::fetchCopcStatisticsEvlrData()
+{
+  uint64_t offset = mLazInfo->firstEvlrOffset();
+  uint32_t evlrCount = mLazInfo->evlrCount();
+
+  QByteArray statisticsEvlrData;
+
+  for ( uint32_t i = 0; i < evlrCount; ++i )
+  {
+    lazperf::evlr_header header;
+    mCopcFile.seekg( offset );
+    char buffer[60];
+    mCopcFile.read( buffer, 60 );
+    header.fill( buffer, 60 );
+
+    // UserID: "qgis", record id: 0
+    if ( header.user_id == "qgis" && header.record_id == 0 )
+    {
+      statisticsEvlrData = QByteArray( header.data_length, Qt::Initialization::Uninitialized );
+      mCopcFile.read( statisticsEvlrData.data(), header.data_length );
+      break;
+    }
+
+    offset += 60 + header.data_length;
+  }
+
+  return statisticsEvlrData;
+}
+
 ///@endcond
diff --git a/src/core/pointcloud/qgscopcpointcloudindex.h b/src/core/pointcloud/qgscopcpointcloudindex.h
index d1f4598d3de..717d86e66d7 100644
--- a/src/core/pointcloud/qgscopcpointcloudindex.h
+++ b/src/core/pointcloud/qgscopcpointcloudindex.h
@@ -32,6 +32,7 @@
 #include "qgspointcloudattribute.h"
 #include "qgsstatisticalsummary.h"
 #include "qgis_sip.h"
+#include "qgspointcloudstatistics.h"
 
 #include "qgslazinfo.h"
 #include "lazperf/vlr.hpp"
@@ -67,6 +68,20 @@ class CORE_EXPORT QgsCopcPointCloudIndex: public QgsPointCloudIndex
     bool isValid() const override;
     QgsPointCloudIndex::AccessType accessType() const override { return QgsPointCloudIndex::Local; };
 
+    /**
+     * Writes the statistics object \a stats into the COPC dataset as an Extended Variable Length Record (EVLR).
+     * Returns true if the data was written successfully.
+     * \since QGIS 3.26
+     */
+    bool writeStatistics( QgsPointCloudStatistics &stats );
+
+    /**
+     * Returns the statistics object contained in the COPC dataset.
+     * If the dataset doesn't contain statistics EVLR, an object with 0 samples will be returned.
+     * \since QGIS 3.26
+     */
+    QgsPointCloudStatistics readStatistics();
+
     /**
      * Copies common properties to the \a destination index
      * \since QGIS 3.26
@@ -86,6 +101,8 @@ class CORE_EXPORT QgsCopcPointCloudIndex: public QgsPointCloudIndex
      */
     virtual void fetchHierarchyPage( uint64_t offset, uint64_t byteSize ) const;
 
+    QByteArray fetchCopcStatisticsEvlrData();
+
     bool mIsValid = false;
     QString mFileName;
     mutable std::ifstream mCopcFile;
diff --git a/src/core/pointcloud/qgslazinfo.h b/src/core/pointcloud/qgslazinfo.h
index 4d5b057a24d..0ea5ea06580 100644
--- a/src/core/pointcloud/qgslazinfo.h
+++ b/src/core/pointcloud/qgslazinfo.h
@@ -102,6 +102,11 @@ class CORE_EXPORT QgsLazInfo
     //! Returns the number of extrabytes contained in the LAZ dataset
     int extrabytesCount() const { return mHeader.ebCount(); }
 
+    //! Returns the absolute offset to the first extended point record in the LAZ file
+    uint64_t firstEvlrOffset() const { return mHeader.evlr_offset; }
+    //! Returns the absolute offset to the first variable length record in the LAZ file
+    uint32_t evlrCount() const { return mHeader.evlr_count; }
+
     //! Returns the coordinate system stored in the LAZ file
     QgsCoordinateReferenceSystem crs() const { return mCrs; }
 
@@ -125,6 +130,11 @@ class CORE_EXPORT QgsLazInfo
     //! Static function to create a QgsLazInfo class from a file over network
     static QgsLazInfo fromUrl( QUrl &url );
 
+#ifndef SIP_RUN
+    //! Returns the LAZPERF header object
+    lazperf::header14 header() const { return mHeader; }
+#endif
+
   private:
     void parseHeader( lazperf::header14 &header );
     void parseCrs();
diff --git a/src/core/pointcloud/qgspointcloudlayer.cpp b/src/core/pointcloud/qgspointcloudlayer.cpp
index d2537e04ee4..e491ccba4fd 100644
--- a/src/core/pointcloud/qgspointcloudlayer.cpp
+++ b/src/core/pointcloud/qgspointcloudlayer.cpp
@@ -39,6 +39,9 @@
 #include "qgsmessagelog.h"
 #include "qgstaskmanager.h"
 #include "qgspointcloudlayerprofilegenerator.h"
+#ifdef HAVE_COPC
+#include "qgscopcpointcloudindex.h"
+#endif
 
 #include <QUrl>
 
@@ -785,6 +788,23 @@ void QgsPointCloudLayer::calculateStatistics()
     QgsMessageLog::logMessage( QObject::tr( "A statistics calculation task for the point cloud %1 is already in progress" ).arg( this->name() ) );
     return;
   }
+#ifdef HAVE_COPC
+  if ( mDataProvider && mDataProvider->index() && mDataProvider->index()->isValid() )
+  {
+    if ( QgsCopcPointCloudIndex *index = qobject_cast<QgsCopcPointCloudIndex *>( mDataProvider->index() ) )
+    {
+      mStatistics = index->readStatistics();
+    }
+  }
+#endif
+  if ( mStatistics.sampledPointsCount() != 0 )
+  {
+    mStatisticsCalculationState = QgsPointCloudLayer::PointCloudStatisticsCalculationState::Calculated;
+    emit statisticsCalculationStateChanged( mStatisticsCalculationState );
+    resetRenderer();
+    return;
+  }
+
   QVector<QgsPointCloudAttribute> attributes = mDataProvider->attributes().attributes();
   // Do not calculate stats for X, Y & Z since the point cloud index contains that
   for ( int i = 0; i < attributes.size(); ++i )
@@ -835,6 +855,15 @@ void QgsPointCloudLayer::calculateStatistics()
     emit statisticsCalculationStateChanged( mStatisticsCalculationState );
     resetRenderer();
     mStatsCalculationTask = 0;
+#ifdef HAVE_COPC
+    if ( mDataProvider && mDataProvider->index() && mDataProvider->index()->isValid() && mDataProvider->name() == QStringLiteral( "pdal" ) && mStatistics.sampledPointsCount() != 0 )
+    {
+      if ( QgsCopcPointCloudIndex *index = qobject_cast<QgsCopcPointCloudIndex *>( mDataProvider->index() ) )
+      {
+        index->writeStatistics( mStatistics );
+      }
+    }
+#endif
   } );
 
   // In case the statistics calculation fails, QgsTask::taskTerminated will be called
@@ -857,7 +886,7 @@ void QgsPointCloudLayer::resetRenderer()
   {
     calculateStatistics();
   }
-  if ( mRenderer->type() == QLatin1String( "extent" ) )
+  if ( !mRenderer || mRenderer->type() == QLatin1String( "extent" ) )
   {
     setRenderer( QgsPointCloudRendererRegistry::defaultRenderer( this ) );
   }
@@ -865,3 +894,5 @@ void QgsPointCloudLayer::resetRenderer()
 
   emit rendererChanged();
 }
+
+
diff --git a/src/core/pointcloud/qgspointcloudstatistics.cpp b/src/core/pointcloud/qgspointcloudstatistics.cpp
index 5e527d55f54..39281f58abc 100644
--- a/src/core/pointcloud/qgspointcloudstatistics.cpp
+++ b/src/core/pointcloud/qgspointcloudstatistics.cpp
@@ -18,8 +18,11 @@
 #include "qgspointcloudstatistics.h"
 
 #include <limits>
+#include <QJsonObject>
+#include <QJsonDocument>
 
 #include "qgspointcloudattribute.h"
+#include "qgsmessagelog.h"
 
 // QgsPointCloudAttributeStatistics
 
@@ -139,3 +142,82 @@ void QgsPointCloudStatistics::combineWith( const QgsPointCloudStatistics &stats
   }
   mSampledPointsCount += stats.mSampledPointsCount;
 }
+
+QByteArray QgsPointCloudStatistics::toStatisticsJson() const
+{
+  QJsonObject obj;
+  obj.insert( QStringLiteral( "sampled-points" ), QJsonValue::fromVariant( sampledPointsCount() ) );
+  QJsonObject stats;
+  for ( const QString &attr : mStatisticsMap.keys() )
+  {
+    QgsPointCloudAttributeStatistics stat = mStatisticsMap.value( attr );
+    stats.insert( attr, attributeStatisticsToJson( stat ) );
+  }
+  obj.insert( QStringLiteral( "stats" ), stats );
+
+  QJsonDocument statsDoc( obj );
+  return statsDoc.toJson( QJsonDocument::Compact );
+}
+
+QgsPointCloudStatistics QgsPointCloudStatistics::fromStatisticsJson( QByteArray statsByteArray )
+{
+  QJsonParseError error;
+  QJsonDocument document = QJsonDocument::fromJson( statsByteArray, &error );
+  if ( error.error != QJsonParseError::NoError )
+  {
+    QgsMessageLog::logMessage( QObject::tr( "Failed to load statistics JSON from COPC file, reason: %1" ).arg( error.errorString() ) );
+    return QgsPointCloudStatistics();
+  }
+
+  QJsonObject statsJson = document.object();
+
+  QgsPointCloudStatistics stats;
+  stats.mSampledPointsCount = statsJson.value( QStringLiteral( "sampled-points" ) ).toInt();
+  if ( statsJson.contains( QStringLiteral( "stats" ) ) )
+  {
+    QJsonObject statsObj = statsJson.value( QStringLiteral( "stats" ) ).toObject();
+    for ( const QString &attr : statsObj.keys() )
+    {
+      QJsonObject obj = statsObj.value( attr ).toObject();
+      QgsPointCloudAttributeStatistics attrStats = fromAttributeStatisticsJson( obj );
+      attrStats.count = stats.mSampledPointsCount;
+      stats.mStatisticsMap.insert( attr, attrStats );
+    }
+  }
+  return stats;
+}
+
+QJsonObject QgsPointCloudStatistics::attributeStatisticsToJson( const QgsPointCloudAttributeStatistics &stats )
+{
+  QJsonObject obj;
+  obj.insert( QStringLiteral( "minimum" ), stats.minimum );
+  obj.insert( QStringLiteral( "maximum" ), stats.maximum );
+  obj.insert( QStringLiteral( "mean" ), stats.mean );
+  if ( !std::isnan( stats.stDev ) )
+  {
+    obj.insert( QStringLiteral( "standard-deviation" ), stats.stDev );
+  }
+  QJsonObject classCount;
+  for ( const int &c : stats.classCount.keys() )
+  {
+    classCount.insert( QString::number( c ), stats.classCount[c] );
+  }
+  obj.insert( QStringLiteral( "class-count" ), classCount );
+  return obj;
+}
+
+QgsPointCloudAttributeStatistics QgsPointCloudStatistics::fromAttributeStatisticsJson( QJsonObject &statsJson )
+{
+  QgsPointCloudAttributeStatistics statsObj;
+  QVariantMap m = statsJson.toVariantMap();
+  statsObj.minimum = m.value( QStringLiteral( "minimum" ), std::numeric_limits<double>::max() ).toDouble();
+  statsObj.maximum = m.value( QStringLiteral( "maximum" ), std::numeric_limits<double>::lowest() ).toDouble();
+  statsObj.mean = m.value( QStringLiteral( "mean" ), 0 ).toDouble();
+  statsObj.stDev = m.value( QStringLiteral( "standard-deviation" ), std::numeric_limits<double>::quiet_NaN() ).toDouble();
+  QJsonObject classCountJson = statsJson.value( QStringLiteral( "class-count" ) ).toObject();
+  for ( const QString &key : classCountJson.keys() )
+  {
+    statsObj.classCount.insert( key.toInt(), classCountJson.value( key ).toInt() );
+  }
+  return statsObj;
+}
diff --git a/src/core/pointcloud/qgspointcloudstatistics.h b/src/core/pointcloud/qgspointcloudstatistics.h
index fd80cfe5b9d..ee2306f233a 100644
--- a/src/core/pointcloud/qgspointcloudstatistics.h
+++ b/src/core/pointcloud/qgspointcloudstatistics.h
@@ -120,6 +120,12 @@ class CORE_EXPORT QgsPointCloudStatistics
     //! Merges the current statistics with the statistics from \a stats
     void combineWith( const QgsPointCloudStatistics &stats );
 
+    //! Converts the current statistics object into JSON object
+    QByteArray toStatisticsJson() const;
+
+    //! Creates a statistics object from the JSON object \a stats
+    static QgsPointCloudStatistics fromStatisticsJson( QByteArray stats );
+
 #ifndef SIP_RUN
     //! Returns a map object containing all the statistics
     QMap<QString, QgsPointCloudAttributeStatistics> statisticsMap() const { return mStatisticsMap; };
@@ -127,6 +133,12 @@ class CORE_EXPORT QgsPointCloudStatistics
   private:
     int mSampledPointsCount = 0;
     QMap<QString, QgsPointCloudAttributeStatistics> mStatisticsMap;
+
+    //! Converts statistics object \a stats into a JSON object
+    static QJsonObject attributeStatisticsToJson( const QgsPointCloudAttributeStatistics &stats );
+
+    //! Creates a statistics object from the JSON object \a stats
+    static QgsPointCloudAttributeStatistics fromAttributeStatisticsJson( QJsonObject &stats );
 };
 
 #endif // QGSPOINTCLOUDSTATISTICS_H
diff --git a/src/core/pointcloud/qgspointcloudstatscalculator.cpp b/src/core/pointcloud/qgspointcloudstatscalculator.cpp
index ceb5dfd7f73..5cefa09d734 100644
--- a/src/core/pointcloud/qgspointcloudstatscalculator.cpp
+++ b/src/core/pointcloud/qgspointcloudstatscalculator.cpp
@@ -175,8 +175,6 @@ QgsPointCloudStatsCalculator::QgsPointCloudStatsCalculator( QgsPointCloudIndex *
 
 }
 
-#include <QThread>
-
 bool QgsPointCloudStatsCalculator::calculateStats( QgsFeedback *feedback, const QVector<QgsPointCloudAttribute> &attributes, qint64 pointsLimit )
 {
   if ( !mIndex->isValid() )
diff --git a/tests/src/providers/testqgscopcprovider.cpp b/tests/src/providers/testqgscopcprovider.cpp
index 6de9d39e957..b4a7926a4f1 100644
--- a/tests/src/providers/testqgscopcprovider.cpp
+++ b/tests/src/providers/testqgscopcprovider.cpp
@@ -46,6 +46,7 @@
 #include "qgsstatisticalsummary.h"
 #include "qgsfeedback.h"
 #include "qgsrangerequestcache.h"
+#include "qgscopcpointcloudindex.h"
 
 /**
  * \ingroup UnitTests
@@ -80,6 +81,7 @@ class TestQgsCopcProvider : public QObject
     void testExtraBytesAttributesValues();
     void testPointCloudIndex();
     void testStatsCalculator();
+    void testSaveLoadStats();
 
     void testQgsRangeRequestCache();
 
@@ -613,7 +615,6 @@ void TestQgsCopcProvider::testExtraBytesAttributesValues()
     {
       for ( const QString &k : keys )
       {
-        qDebug() << i << k << identifiedPoints[i][k].toDouble() << " " << expectedPoints[i][k].toDouble();
         QCOMPARE( identifiedPoints[i][k].toDouble(), expectedPoints[i][k].toDouble() );
       }
     }
@@ -891,5 +892,34 @@ void TestQgsCopcProvider::testQgsRangeRequestCache()
   }
 }
 
+void TestQgsCopcProvider::testSaveLoadStats()
+{
+  QgsPointCloudStatistics calculatedStats;
+  QgsPointCloudStatistics readStats;
+  {
+    std::unique_ptr< QgsPointCloudLayer > layer = std::make_unique< QgsPointCloudLayer >( mTestDataDir + QStringLiteral( "point_clouds/copc/lone-star.copc.laz" ), QStringLiteral( "layer" ), QStringLiteral( "copc" ) );
+    QVERIFY( layer->isValid() );
+
+    QVERIFY( layer->dataProvider() && layer->dataProvider()->isValid() && layer->dataProvider()->index() );
+    QgsCopcPointCloudIndex *index = qobject_cast<QgsCopcPointCloudIndex *>( layer->dataProvider()->index() );
+
+    calculatedStats = layer->statistics();
+    index->writeStatistics( calculatedStats );
+  }
+
+  {
+    std::unique_ptr< QgsPointCloudLayer > layer = std::make_unique< QgsPointCloudLayer >( mTestDataDir + QStringLiteral( "point_clouds/copc/lone-star.copc.laz" ), QStringLiteral( "layer" ), QStringLiteral( "copc" ) );
+    QVERIFY( layer->isValid() );
+
+    QVERIFY( layer->dataProvider() && layer->dataProvider()->isValid() && layer->dataProvider()->index() );
+
+    QgsCopcPointCloudIndex *index = qobject_cast<QgsCopcPointCloudIndex *>( layer->dataProvider()->index() );
+    readStats = index->readStatistics();
+  }
+
+  QVERIFY( calculatedStats.sampledPointsCount() == readStats.sampledPointsCount() );
+  QVERIFY( calculatedStats.toStatisticsJson() == readStats.toStatisticsJson() );
+}
+
 QGSTEST_MAIN( TestQgsCopcProvider )
 #include "testqgscopcprovider.moc"