Branch data     Line data    Source code

       1                 :            : /***************************************************************************
       2                 :            :     qgsclassificationjenks.h
       3                 :            :     ---------------------
       4                 :            :     begin                : September 2019
       5                 :            :     copyright            : (C) 2019 by Denis Rouzaud
       6                 :            :     email                : denis@opengis.ch
       7                 :            :  ***************************************************************************
       8                 :            :  *                                                                         *
       9                 :            :  *   This program is free software; you can redistribute it and/or modify  *
      10                 :            :  *   it under the terms of the GNU General Public License as published by  *
      11                 :            :  *   the Free Software Foundation; either version 2 of the License, or     *
      12                 :            :  *   (at your option) any later version.                                   *
      13                 :            :  *                                                                         *
      14                 :            :  ***************************************************************************/
      15                 :            : 
      16                 :            : #include <limits>
      17                 :            : #include "qgsclassificationjenks.h"
      18                 :            : #include "qgsapplication.h"
      19                 :            : 
      20                 :            : #if QT_VERSION >= QT_VERSION_CHECK(5, 15, 0)
      21                 :            : #include <QRandomGenerator>
      22                 :            : #endif
      23                 :            : 
      24                 :          5 : QgsClassificationJenks::QgsClassificationJenks()
      25                 :          5 :   : QgsClassificationMethod()
      26                 :         10 : {
      27                 :            : 
      28                 :          5 : }
      29                 :            : 
      30                 :          0 : QString QgsClassificationJenks::name() const
      31                 :            : {
      32                 :          0 :   return QObject::tr( "Natural Breaks (Jenks)" );
      33                 :            : }
      34                 :            : 
      35                 :         10 : QString QgsClassificationJenks::id() const
      36                 :            : {
      37                 :         20 :   return QStringLiteral( "Jenks" );
      38                 :            : }
      39                 :            : 
      40                 :          0 : QgsClassificationMethod *QgsClassificationJenks::clone() const
      41                 :            : {
      42                 :          5 :   QgsClassificationJenks *c = new QgsClassificationJenks();
      43                 :          0 :   copyBase( c );
      44                 :          0 :   return c;
      45                 :          0 : }
      46                 :            : 
      47                 :          0 : QIcon QgsClassificationJenks::icon() const
      48                 :            : {
      49                 :          0 :   return QgsApplication::getThemeIcon( "classification_methods/mClassificationNaturalBreak.svg" );
      50                 :          0 : }
      51                 :            : 
      52                 :            : 
      53                 :          0 : QList<double> QgsClassificationJenks::calculateBreaks( double &minimum, double &maximum,
      54                 :            :     const QList<double> &values, int nclasses )
      55                 :            : {
      56                 :            :   // Jenks Optimal (Natural Breaks) algorithm
      57                 :            :   // Based on the Jenks algorithm from the 'classInt' package available for
      58                 :            :   // the R statistical prgramming language, and from Python code from here:
      59                 :            :   // http://danieljlewis.org/2010/06/07/jenks-natural-breaks-algorithm-in-python/
      60                 :            :   // and is based on a JAVA and Fortran code available here:
      61                 :            :   // https://stat.ethz.ch/pipermail/r-sig-geo/2006-March/000811.html
      62                 :            : 
      63                 :            :   // Returns class breaks such that classes are internally homogeneous while
      64                 :            :   // assuring heterogeneity among classes.
      65                 :            : 
      66                 :          0 :   if ( values.isEmpty() )
      67                 :          0 :     return QList<double>();
      68                 :            : 
      69                 :          0 :   if ( nclasses <= 1 )
      70                 :            :   {
      71                 :          0 :     return QList<double>() <<  maximum;
      72                 :            :   }
      73                 :            : 
      74                 :          0 :   if ( nclasses >= values.size() )
      75                 :            :   {
      76                 :          0 :     return values;
      77                 :            :   }
      78                 :            : 
      79                 :          0 :   QVector<double> sample;
      80                 :            : 
      81                 :            :   // if we have lots of values, we need to take a random sample
      82                 :          0 :   if ( values.size() > mMaximumSize )
      83                 :            :   {
      84                 :            :     // for now, sample at least maximumSize values or a 10% sample, whichever
      85                 :            :     // is larger. This will produce a more representative sample for very large
      86                 :            :     // layers, but could end up being computationally intensive...
      87                 :            : 
      88                 :          0 :     sample.resize( std::max( mMaximumSize, static_cast<int>( values.size() ) / 10 ) );
      89                 :            : 
      90                 :          0 :     QgsDebugMsgLevel( QStringLiteral( "natural breaks (jenks) sample size: %1" ).arg( sample.size() ), 2 );
      91                 :          0 :     QgsDebugMsgLevel( QStringLiteral( "values:%1" ).arg( values.size() ), 2 );
      92                 :            : 
      93                 :          0 :     sample[ 0 ] = minimum;
      94                 :          0 :     sample[ 1 ] = maximum;
      95                 :            : 
      96                 :          0 :     for ( int i = 2; i < sample.size(); i++ )
      97                 :            :     {
      98                 :            :       // pick a random integer from 0 to n
      99                 :            : 
     100                 :            : #if QT_VERSION >= QT_VERSION_CHECK(5, 15, 0)
     101                 :          0 :       double r = QRandomGenerator::global()->generate();
     102                 :            : #else
     103                 :            :       double r = qrand();
     104                 :            : #endif
     105                 :          0 :       int j = std::floor( r / RAND_MAX * ( values.size() - 1 ) );
     106                 :          0 :       sample[ i ] = values[ j ];
     107                 :          0 :     }
     108                 :          0 :   }
     109                 :            :   else
     110                 :            :   {
     111                 :          0 :     sample = values.toVector();
     112                 :            :   }
     113                 :            : 
     114                 :          0 :   int n = sample.size();
     115                 :            : 
     116                 :            :   // sort the sample values
     117                 :          0 :   std::sort( sample.begin(), sample.end() );
     118                 :            : 
     119                 :          0 :   QVector< QVector<int> > matrixOne( n + 1 );
     120                 :          0 :   QVector< QVector<double> > matrixTwo( n + 1 );
     121                 :            : 
     122                 :          0 :   for ( int i = 0; i <= n; i++ )
     123                 :            :   {
     124                 :          0 :     matrixOne[i].resize( nclasses + 1 );
     125                 :          0 :     matrixTwo[i].resize( nclasses + 1 );
     126                 :          0 :   }
     127                 :            : 
     128                 :          0 :   for ( int i = 1; i <= nclasses; i++ )
     129                 :            :   {
     130                 :          0 :     matrixOne[0][i] = 1;
     131                 :          0 :     matrixOne[1][i] = 1;
     132                 :          0 :     matrixTwo[0][i] = 0.0;
     133                 :          0 :     for ( int j = 2; j <= n; j++ )
     134                 :            :     {
     135                 :          0 :       matrixTwo[j][i] = std::numeric_limits<double>::max();
     136                 :          0 :     }
     137                 :          0 :   }
     138                 :            : 
     139                 :          0 :   for ( int l = 2; l <= n; l++ )
     140                 :            :   {
     141                 :          0 :     double s1 = 0.0;
     142                 :          0 :     double s2 = 0.0;
     143                 :          0 :     int w = 0;
     144                 :            : 
     145                 :          0 :     double v = 0.0;
     146                 :            : 
     147                 :          0 :     for ( int m = 1; m <= l; m++ )
     148                 :            :     {
     149                 :          0 :       int i3 = l - m + 1;
     150                 :            : 
     151                 :          0 :       double val = sample[ i3 - 1 ];
     152                 :            : 
     153                 :          0 :       s2 += val * val;
     154                 :          0 :       s1 += val;
     155                 :          0 :       w++;
     156                 :            : 
     157                 :          0 :       v = s2 - ( s1 * s1 ) / static_cast< double >( w );
     158                 :          0 :       int i4 = i3 - 1;
     159                 :          0 :       if ( i4 != 0 )
     160                 :            :       {
     161                 :          0 :         for ( int j = 2; j <= nclasses; j++ )
     162                 :            :         {
     163                 :          0 :           if ( matrixTwo[l][j] >= v + matrixTwo[i4][j - 1] )
     164                 :            :           {
     165                 :          0 :             matrixOne[l][j] = i4;
     166                 :          0 :             matrixTwo[l][j] = v + matrixTwo[i4][j - 1];
     167                 :          0 :           }
     168                 :          0 :         }
     169                 :          0 :       }
     170                 :          0 :     }
     171                 :          0 :     matrixOne[l][1] = 1;
     172                 :          0 :     matrixTwo[l][1] = v;
     173                 :          0 :   }
     174                 :            : 
     175                 :          0 :   QVector<double> breaks( nclasses );
     176                 :          0 :   breaks[nclasses - 1] = sample[n - 1];
     177                 :            : 
     178                 :          0 :   for ( int j = nclasses, k = n; j >= 2; j-- )
     179                 :            :   {
     180                 :          0 :     int id = matrixOne[k][j] - 1;
     181                 :          0 :     breaks[j - 2] = sample[id];
     182                 :          0 :     k = matrixOne[k][j] - 1;
     183                 :          0 :   }
     184                 :            : 
     185                 :          0 :   return breaks.toList();
     186                 :          0 : }