<div dir="ltr">Hello,<div>Sorry if this gets posted twice I previously attached a .gpkg and it bounced due to the 100kb limit. <br><div>I am enjoying this fascinating discussion. This is all fairly new to me so I looked some things up and want to make sure that I am getting things right, </div><div><ol><li style="margin-left:15px">In QGIS the $area function uses the sphere/ellipsoid settings of the project not of the CRS of the layer. </li><li style="margin-left:15px">QGIS uses the PROJ library to calculate area</li><li style="margin-left:15px">According to the PROJ documentation <a href="http://even.rouault.free.fr/proj_cpp_api/rst_generated/html/geodesic.html" target="_blank">even.rouault.free.fr/proj_cpp_api/rst_generated/html/geodesic.html</a> PROJ uses GeographicLib to calculate area. </li></ol><div>If this is all true then if as  Jukka Rahkonen describes the layer and the project were both EPSG:4326 and the ellipsoid in the Project-Settings was WGS 84 (EPSG:7030) then the QGIS measurements should be consistent (if they are correct is open to discussion).with this GeographicLib web calculator that says it uses the WGS84 ellipsoid <a href="http://geographiclib.sourceforge.io/cgi-bin/Planimeter?type=polygon&rhumb=geodesic&input=59.95688345+20.13293641%0D%0A60.47397663+26.94617837%0D%0A62.56499443+29.74782155%0D%0A68.70650340+27.45254202%0D%0A68.24937206+23.75771765%0D%0A65.27444593+25.42698984%0D%0A63.10353609+21.51545237%0D%0A61.12318104+21.40562760%0D%0A60.40477513+19.41123592%0D%0A59.95688345+20.13293641&option=Submit" target="_blank">geographiclib.sourceforge.io/cgi-bin/Planimeter?type=polygon&rhumb=geodesic&input=59.95688345+20.13293641%0D%0A60.47397663+26.94617837%0D%0A62.56499443+29.74782155%0D%0A68.70650340+27.45254202%0D%0A68.24937206+23.75771765%0D%0A65.27444593+25.42698984%0D%0A63.10353609+21.51545237%0D%0A61.12318104+21.40562760%0D%0A60.40477513+19.41123592%0D%0A59.95688345+20.13293641&option=Submit</a><br></div></div><div><br></div><div>Using Jukka Rahkonen's points and processes I ended up with the same non consistent results as they reported. <strike>Attached is a geopackage with the polygon that I used,</strike></div></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Wed, Jan 6, 2021 at 10:29 AM Sebastian Gutwein <<a href="mailto:bas@rdgland.com">bas@rdgland.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div dir="ltr">Hello,<div>I am enjoying this fascinating discussion. This is all fairly new to me so I looked some things up and want to make sure that I am getting things right, </div><div><ol><li>In QGIS the $area function uses the sphere/ellipsoid settings of the project not of the CRS of the layer. </li><li>QGIS uses the PROJ library to calculate area</li><li>According to the PROJ documentation <a href="http://even.rouault.free.fr/proj_cpp_api/rst_generated/html/geodesic.html" target="_blank">even.rouault.free.fr/proj_cpp_api/rst_generated/html/geodesic.html</a> PROJ uses GeographicLib to calculate area. </li></ol><div>If this is all true then if as 

Jukka Rahkonen describes the layer and the project were both EPSG:4326 and the
              ellipsoid in the Project-Settings was WGS 84 (EPSG:7030) then the QGIS measurements should be consistent (if they are correct is open to discussion).with this GeographicLib web calculator that says it uses the WGS84 ellipsoid <a href="http://geographiclib.sourceforge.io/cgi-bin/Planimeter?type=polygon&rhumb=geodesic&input=59.95688345+20.13293641%0D%0A60.47397663+26.94617837%0D%0A62.56499443+29.74782155%0D%0A68.70650340+27.45254202%0D%0A68.24937206+23.75771765%0D%0A65.27444593+25.42698984%0D%0A63.10353609+21.51545237%0D%0A61.12318104+21.40562760%0D%0A60.40477513+19.41123592%0D%0A59.95688345+20.13293641&option=Submit" target="_blank">geographiclib.sourceforge.io/cgi-bin/Planimeter?type=polygon&rhumb=geodesic&input=59.95688345+20.13293641%0D%0A60.47397663+26.94617837%0D%0A62.56499443+29.74782155%0D%0A68.70650340+27.45254202%0D%0A68.24937206+23.75771765%0D%0A65.27444593+25.42698984%0D%0A63.10353609+21.51545237%0D%0A61.12318104+21.40562760%0D%0A60.40477513+19.41123592%0D%0A59.95688345+20.13293641&option=Submit</a><br></div></div><div><br></div><div>Using Jukka Rahkonen's points and processes I ended up with the same non consistent results as they reported. Attached is a geopackage with the polygon that I used, </div><div>-Sebastian</div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Tue, Jan 5, 2021 at 3:52 PM Nicolas Cadieux <<a href="mailto:njacadieux.gitlab@gmail.com" target="_blank">njacadieux.gitlab@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
  
    
  
  <div>
    <p>Hi,</p>
    <p>QGIS is currently built on Proj version 6.3.2-1.  If the other
      libraries are using a spheroid by default, then they use a sphere
      for speed and not an ellipsoid for precision.  You can probably
      force this measurement in QGIS by creating a custom CRS with a
      spheroid rather than an ellipsoid. The other option is to use the
      python in QGIS and to force a geoid.<br>
    </p>
    <p>Nicolas   <br>
    </p>
    <div>On 2021-01-05 1:33 p.m., Rahkonen Jukka
      (MML) wrote:<br>
    </div>
    <blockquote type="cite">
      
      
      
      <div>
        <p class="MsoNormal">Hi,<u></u><u></u></p>
        <p class="MsoNormal"><u></u> <u></u></p>
        <p class="MsoNormal"><span lang="EN-US">I suppose that PostGIS
            is using the WGS 84 ellipsoid but I am not sure where from
            the documentation I could find that information. The ST_Area
            document
            <a href="https://postgis.net/docs/ST_Area.html" target="_blank">https://postgis.net/docs/ST_Area.html</a>
            says only “For geography types by default area is determined
            on a spheroid with units in square meters”. Same thing with
            Spatialite, documentation suggests just that it is “the”
            spheroid
            <a href="http://www.gaia-gis.it/gaia-sins/spatialite-sql-latest.html" target="_blank">http://www.gaia-gis.it/gaia-sins/spatialite-sql-latest.html</a>.<u></u><u></u></span></p>
        <p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p>
        <p class="MsoNormal"><span lang="EN-US">I did not notice this
            paragraph in the ST_Area document earlier “Enhanced: 2.2.0 -
            measurement on spheroid performed with GeographicLib for
            improved accuracy and robustness. Requires Proj >= 4.9.0
            to take advantage of the new feature.” So no wonder that the
            web app and PostGIS give the same results because they both
            use GeographicLib. And SpatiaLite 5.0 is using RTTopo that
            is a library that is based on LWGeom so close connection in
            there too.<u></u><u></u></span></p>
        <p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p>
        <p class="MsoNormal"><span lang="EN-US">So perhaps the way to
            get identical areas from QGIS would be to make it to use
            GeographicLib as well. I have no idea if it is a realistic
            approach and worth making a feature request.<u></u><u></u></span></p>
        <p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p>
        <p class="MsoNormal"><span lang="EN-US">-Jukka Rahkonen-<u></u><u></u></span></p>
        <p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p>
        <div>
          <div style="border-right:none;border-bottom:none;border-left:none;border-top:1pt solid rgb(225,225,225);padding:3pt 0cm 0cm">
            <p class="MsoNormal"><b><span>Lähettäjä:</span></b><span> Nicolas Cadieux
                <a href="mailto:njacadieux.gitlab@gmail.com" target="_blank"><njacadieux.gitlab@gmail.com></a>
                <br>
                <b>Lähetetty:</b> tiistai 5. tammikuuta 2021 20.04<br>
                <b>Vastaanottaja:</b> Rahkonen Jukka (MML)
                <a href="mailto:jukka.rahkonen@maanmittauslaitos.fi" target="_blank"><jukka.rahkonen@maanmittauslaitos.fi></a>;
                <a href="mailto:qgis-user@lists.osgeo.org" target="_blank">qgis-user@lists.osgeo.org</a><br>
                <b>Aihe:</b> Re: [Qgis-user] Confusion with ellipsoidal
                method of $area<u></u><u></u></span></p>
          </div>
        </div>
        <p class="MsoNormal"><u></u> <u></u></p>
        <p>Hi,<span><u></u><u></u></span></p>
        <p>Your method in QGIS is sound.  Area is calculated using the
          wgs84 ellipsoid EPSG 7030. If you reproject to wgs84 zone 35N
          (I think this is close), area goes from 249566957499.7546m2 
          to 249566957499.721m2 or a difference of 0.0336 m2.  I don't
          think densification would change things much.<u></u><u></u></p>
        <p>My question is the following: You know that QGIS uses the
          WGS84 Ellipsoid.  What Ellipsoid are using used in the other
          software???  If you don't know, I would proceed until you
          figure that out.<u></u><u></u></p>
        <p>This is the currently used ellipsoid in the proj database in
          QGIS 3.16.  You can get this by typing proj -le in the OSFeo4W
          Shell.<u></u><u></u></p>
        <p><br>
              MERIT a=6378137.0      rf=298.257       MERIT 1983<br>
              SGS85 a=6378136.0      rf=298.257       Soviet Geodetic
          System 85<br>
              GRS80 a=6378137.0      rf=298.257222101 GRS 1980(IUGG,
          1980)<br>
              IAU76 a=6378140.0      rf=298.257       IAU 1976<br>
               airy a=6377563.396    rf=299.3249646   Airy 1830<br>
             APL4.9 a=6378137.0      rf=298.25        Appl. Physics.
          1965<br>
              NWL9D a=6378145.0      rf=298.25        Naval Weapons
          Lab., 1965<br>
           mod_airy a=6377340.189    b=6356034.446    Modified Airy<br>
             andrae a=6377104.43     rf=300.0         Andrae 1876 (Den.,
          Iclnd.)<br>
             danish a=6377019.2563   rf=300.0         Andrae 1876
          (Denmark, Iceland)<br>
            aust_SA a=6378160.0      rf=298.25        Australian Natl
          & S. Amer. 1969<br>
              GRS67 a=6378160.0      rf=298.2471674270 GRS 67(IUGG 1967)<br>
            GSK2011 a=6378136.5      rf=298.2564151   GSK-2011<br>
             bessel a=6377397.155    rf=299.1528128   Bessel 1841<br>
           bess_nam a=6377483.865    rf=299.1528128   Bessel 1841
          (Namibia)<br>
             clrk66 a=6378206.4      b=6356583.8      Clarke 1866<br>
             clrk80 a=6378249.145    rf=293.4663      Clarke 1880 mod.<br>
          clrk80ign a=6378249.2      rf=293.4660212936269 Clarke 1880
          (IGN).<br>
                CPM a=6375738.7      rf=334.29        Comm. des Poids et
          Mesures 1799<br>
             delmbr a=6376428.       rf=311.5         Delambre 1810
          (Belgium)<br>
            engelis a=6378136.05     rf=298.2566      Engelis 1985<br>
            evrst30 a=6377276.345    rf=300.8017      Everest 1830<br>
            evrst48 a=6377304.063    rf=300.8017      Everest 1948<br>
            evrst56 a=6377301.243    rf=300.8017      Everest 1956<br>
            evrst69 a=6377295.664    rf=300.8017      Everest 1969<br>
            evrstSS a=6377298.556    rf=300.8017      Everest (Sabah
          & Sarawak)<br>
            fschr60 a=6378166.       rf=298.3         Fischer (Mercury
          Datum) 1960<br>
           fschr60m a=6378155.       rf=298.3         Modified Fischer
          1960<br>
            fschr68 a=6378150.       rf=298.3         Fischer 1968<br>
            helmert a=6378200.       rf=298.3         Helmert 1906<br>
              hough a=6378270.0      rf=297.          Hough<br>
               intl a=6378388.0      rf=297.          International 1909
          (Hayford)<br>
              krass a=6378245.0      rf=298.3         Krassovsky, 1942<br>
              kaula a=6378163.       rf=298.24        Kaula 1961<br>
              lerch a=6378139.       rf=298.257       Lerch 1979<br>
              mprts a=6397300.       rf=191.          Maupertius 1738<br>
           new_intl a=6378157.5      b=6356772.2      New International
          1967<br>
            plessis a=6376523.       b=6355863.       Plessis 1817
          (France)<br>
               PZ90 a=6378136.0      rf=298.25784     PZ-90<br>
             SEasia a=6378155.0      b=6356773.3205   Southeast Asia<br>
            walbeck a=6376896.0      b=6355834.8467   Walbeck<br>
              WGS60 a=6378165.0      rf=298.3         WGS 60<br>
              WGS66 a=6378145.0      rf=298.25        WGS 66<br>
              WGS72 a=6378135.0      rf=298.26        WGS 72<br>
              WGS84 a=6378137.0      rf=298.257223563 WGS 84<br>
             sphere a=6370997.0      b=6370997.0      Normal Sphere
          (r=6370997)<u></u><u></u></p>
        <p>Nicolas<u></u><u></u></p>
        <p><u></u> <u></u></p>
        <div>
          <p class="MsoNormal">On 2021-01-05 11:43 a.m., Rahkonen Jukka
            (MML) wrote:<u></u><u></u></p>
        </div>
        <blockquote style="margin-top:5pt;margin-bottom:5pt">
          <p class="MsoNormal">Hi,<u></u><u></u></p>
          <p class="MsoNormal"> <u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">I wonder what method
              QGIS is using when it computes ellipsoidal area with $area
              function. I made a test with this EPSG:4326 polygon:</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">POLYGON (( 20.13293641
              59.95688345, 26.94617837 60.47397663, 29.74782155
              62.56499443, 27.45254202 68.7065034, 23.75771765
              68.24937206, 25.42698984 65.27444593, 21.51545237
              63.10353609, 21.4056276 61.12318104, 19.41123592
              60.40477513, 20.13293641 59.95688345 ))</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US"> </span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">I checked that the CRS
              of the project and layer were both EPSG:4326. The
              ellipsoid in the Project-Settings was WGS 84 (EPSG:7030).
              The area that $area returns is</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">249566957499.7546</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US"> </span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">As a references I used
              PostGIS and web site
              <a href="https://geographiclib.sourceforge.io/scripts/geod-calc.html#area" target="_blank">https://geographiclib.sourceforge.io/scripts/geod-calc.html#area</a>.
              With my test polygon they both give this result:</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">251199344354.4308</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US"> </span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">PostGIS returns bigger
              area. The difference is 0.654% so not huge but not
              negligible either.</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US"> </span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">A third test with
              SpatiaLite 5.0 gives a result that is very close to
              PostGIS (difference 0.001%).</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">251195856999.549927</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US"> </span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">As a conclusion I have
              decided to trust in PostGIS because it gives the same
              results than the web site
              <a href="https://geographiclib.sourceforge.io/scripts/geod-calc.html" target="_blank">https://geographiclib.sourceforge.io/scripts/geod-calc.html</a>
              that feels scientifically sound. However, I wonder if I
              have used QGIS in a correct way or if there is anything I
              could do for getting areas to match better with PostGIS
              for example by densifying the long segments in my test
              polygon.</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US"> </span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">For getting testing
              easy I attach directly runnable  SQL for PostGIS and
              Spatialite (version 5.0 with RTTopo is needed) as well as
              a coordinate list for the web app.</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US"> </span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">PostGIS</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">select
              st_area(st_geogfromtext('POLYGON ((</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">20.13293641  
              59.95688345,</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">26.94617837  
              60.47397663,</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">29.74782155  
              62.56499443,</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">27.45254202  
              68.70650340,</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">23.75771765  
              68.24937206,</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">25.42698984  
              65.27444593,</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">21.51545237  
              63.10353609,</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">21.40562760  
              61.12318104,</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">19.41123592  
              60.40477513,</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">20.13293641  
              59.95688345))') , true);</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US"> </span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">Spatialite</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">select
              st_area(st_geomfromtext('POLYGON ((</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">20.13293641  
              59.95688345,</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">26.94617837  
              60.47397663,</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">29.74782155  
              62.56499443,</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">27.45254202  
              68.70650340,</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">23.75771765  
              68.24937206,</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">25.42698984  
              65.27444593,</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">21.51545237  
              63.10353609,</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">21.40562760  
              61.12318104,</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">19.41123592  
              60.40477513,</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">20.13293641  
              59.95688345))') ,true);</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US"> </span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">Web site (notice
              lat-lon order)</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">59.95688345
              20.13293641 </span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">60.47397663
              26.94617837 </span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">62.56499443
              29.74782155 </span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">68.70650340
              27.45254202 </span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">68.24937206
              23.75771765 </span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">65.27444593
              25.42698984 </span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">63.10353609
              21.51545237 </span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">61.12318104
              21.40562760 </span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">60.40477513
              19.41123592 </span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">59.95688345
              20.13293641</span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US"> </span><u></u><u></u></p>
          <p class="MsoNormal"><span lang="EN-US">-Jukka Rahkonen- </span><u></u><u></u></p>
          <p class="MsoNormal"><span><br>
              <br>
              <u></u><u></u></span></p>
          <pre>_______________________________________________<u></u><u></u></pre>
          <pre>Qgis-user mailing list<u></u><u></u></pre>
          <pre><a href="mailto:Qgis-user@lists.osgeo.org" target="_blank">Qgis-user@lists.osgeo.org</a><u></u><u></u></pre>
          <pre>List info: <a href="https://lists.osgeo.org/mailman/listinfo/qgis-user" target="_blank">https://lists.osgeo.org/mailman/listinfo/qgis-user</a><u></u><u></u></pre>
          <pre>Unsubscribe: <a href="https://lists.osgeo.org/mailman/listinfo/qgis-user" target="_blank">https://lists.osgeo.org/mailman/listinfo/qgis-user</a><u></u><u></u></pre>
        </blockquote>
        <pre>-- <u></u><u></u></pre>
        <pre>Nicolas Cadieux<u></u><u></u></pre>
        <pre><a href="https://gitlab.com/njacadieux" target="_blank">https://gitlab.com/njacadieux</a><u></u><u></u></pre>
      </div>
    </blockquote>
    <pre cols="72">-- 
Nicolas Cadieux
<a href="https://gitlab.com/njacadieux" target="_blank">https://gitlab.com/njacadieux</a></pre>
  </div>

_______________________________________________<br>
Qgis-user mailing list<br>
<a href="mailto:Qgis-user@lists.osgeo.org" target="_blank">Qgis-user@lists.osgeo.org</a><br>
List info: <a href="https://lists.osgeo.org/mailman/listinfo/qgis-user" rel="noreferrer" target="_blank">https://lists.osgeo.org/mailman/listinfo/qgis-user</a><br>
Unsubscribe: <a href="https://lists.osgeo.org/mailman/listinfo/qgis-user" rel="noreferrer" target="_blank">https://lists.osgeo.org/mailman/listinfo/qgis-user</a><br>
</blockquote></div>
</div>
</blockquote></div>