[R-sig-Geo] Difference between area calculations by QGIS and sf package

[Roger Bivand]

On Mon, 27 Sep 2021, Rik Ferreira wrote:

> Dear Roger, thanks for the explanation.
>
> The error between QGIS and sf drops a lot with sf_use_s2(FALSE) but it's
> still bigger than 0. I believe it is due to different parameters passed
> because GEOS and PROJ are QGIS dependencies.

No, the differences have nothing to do with GEOS, which only handles 
planar geometries, or PROJ, which expresses the coordinate reference 
system (here "EPSG:4674", a GEOGCRS in decimal degrees). The differences 
are coming from difference ways of estimating area on the sphere (s2) or a 
chosen ellipsoid (lwgeom, QGIS) all using different approaches. Try 
projecting to UTM zone 23; then GEOS would be in play, because the 
geometry is planar.

Hope this clarifies,

Roger


>
> I will run more tests to find out what may be causing this issue but the
> problem is solved for now.
>
> Thank you!
>
> Em seg., 27 de set. de 2021 às 12:32, Roger Bivand <Roger.Bivand using nhh.no>
> escreveu:
>
>> On Mon, 27 Sep 2021, Rik Ferreira wrote:
>>
>>> Dear list members,
>>>
>>> Why does sf::st_area returns a different area than the output of QGIS
>> $area
>>> for the same geometry?
>>>
>>> A reproducible example (state of Minas Gerais, Brazil):
>>>
>>> library(sf)
>>>> mg <- geobr::read_state(31)
>>>> st_area(mg)
>>>>
>>>
>>> sf::st_area: 588358849314 [m^2]
>>> QGIS $area: 586521037474,860 m²
>>>
>>> The layer used in QGIS was exported with st_write() in GeoJSON format.
>>>
>>> The error is: 1837811840 m².
>>
>>> st_is_longlat(mg)
>> [1] TRUE
>>> sf_use_s2()
>> [1] TRUE
>>> sf_use_s2(FALSE)
>> Spherical geometry (s2) switched off
>>> (o <- st_area(mg))
>> 586520883040 [m^2]
>>> sf_use_s2(TRUE)
>> Spherical geometry (s2) switched on
>>> (p <- st_area(mg))
>> 588358849314 [m^2]
>>> p - o
>> 1837966275 [m^2]
>>
>> See https://r-spatial.org/r/2020/06/17/s2.html
>>
>> When sf_use_s2() is FALSE, see ?st_area ("ellipsoidal distances are
>> computed using st_geod_distance which uses function 'geod_inverse' from
>> GeographicLib (part of PROJ)"). You'd have to establish what code QGIS
>> uses to know why it and GeographicLib differ.
>>
>> Hope this clarifies
>>
>> Roger
>>
>>
>>>
>>> Thank you for the attention.
>>>
>>> Att.
>>>
>>>
>>
>> --
>> Roger Bivand
>> Emeritus Professor
>> Department of Economics, Norwegian School of Economics,
>> Postboks 3490 Ytre Sandviken, 5045 Bergen, Norway.
>> e-mail: Roger.Bivand using nhh.no
>> https://orcid.org/0000-0003-2392-6140
>> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>
>
>
>

-- 
Roger Bivand
Emeritus Professor
Department of Economics, Norwegian School of Economics,
Postboks 3490 Ytre Sandviken, 5045 Bergen, Norway.
e-mail: Roger.Bivand using nhh.no
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en