[GRASS-user] v.to.db area calculations

GRASS GRASS user
,
1

Hello,

I’m trying to calculate the area of a polygon using v.to.db as follows (through R):

execGRASS(“v.to.db”, flags = c(“p”,“quiet”), parameters = list(map = “test”, units = “meters”, option = “area”, separator = “comma”),intern = TRUE)

I assume that the output of this command will be a meters squared calculation of the area inside the polygon ‘test’.

However, in QGIS, I manually calculated the area the they do not match. QGIS (which I believe to be true) shows 0.128m2, and the output from v.to.db shows 22.42226m2.

I am working in WGS84.

Can anyone shed some light on how v.to.db calculates area and why it might be going wrong please?

Thank you

James

···

James Duffy
PhD Researcher
Environment and Sustainability Institute
Penryn Campus
University of Exeter
Penryn
Cornwall
TR10 9FE

2

Hi James,

QGIS uses a ellipsoidal area measurement, while GRASS (like most of the other GIS) uses a Cartesian measurement (others, please correct me if that is not the case in WGS84).

See amongst others:
https://issues.qgis.org/issues/12057 and
https://issues.qgis.org/issues/14675

Would you mind comparing area measurement of a third software (SAGA, PostGIS, ...)?

Cheers
Stefan

________________________________________
Von: grass-user [grass-user-bounces@lists.osgeo.org] im Auftrag von James Duffy [james.philip.duffy@gmail.com]
Gesendet: Donnerstag, 18. Mai 2017 16:22
An: grass-user
Betreff: [GRASS-user] v.to.db area calculations

Hello,

I'm trying to calculate the area of a polygon using v.to.db as follows (through R):

execGRASS("v.to.db", flags = c("p","quiet"), parameters = list(map = "test", units = "meters", option = "area", separator = "comma"),intern = TRUE)

I assume that the output of this command will be a meters squared calculation of the area inside the polygon 'test'.

However, in QGIS, I manually calculated the area the they do not match. QGIS (which I believe to be true) shows 0.128m2, and the output from v.to.db shows 22.42226m2.

I am working in WGS84.

Can anyone shed some light on how v.to.db calculates area and why it might be going wrong please?

Thank you

James

--
James Duffy
PhD Researcher
Environment and Sustainability Institute
Penryn Campus
University of Exeter
Penryn
Cornwall
TR10 9FE

3

On 18/05/17 18:41, Blumentrath, Stefan wrote:

Hi James,

QGIS uses a ellipsoidal area measurement, while GRASS (like most of the other GIS) uses a Cartesian measurement (others, please correct me if that is not the case in WGS84).

See amongst others:
https://issues.qgis.org/issues/12057 and
https://issues.qgis.org/issues/14675

Would you mind comparing area measurement of a third software (SAGA, PostGIS, ...)?

I would suggest making tests with polygons of which you know the area (e.g. simple squares to begin with)

Moritz

4

Hello,

Thanks for the reply. I don’t have SAGA or PostGIS set up at the moment. (have sent Stefan the shapefile off-list).

Moritz - I know that the QGIS value is true. The polygon has been drawn over a feature on a raster image (which I know is correct in its dimensions).

Thanks

James

···

On 19 May 2017 at 08:36, Moritz Lennert <mlennert@club.worldonline.be> wrote:

On 18/05/17 18:41, Blumentrath, Stefan wrote:

Hi James,

QGIS uses a ellipsoidal area measurement, while GRASS (like most of the other GIS) uses a Cartesian measurement (others, please correct me if that is not the case in WGS84).

See amongst others:
https://issues.qgis.org/issues/12057 and
https://issues.qgis.org/issues/14675

Would you mind comparing area measurement of a third software (SAGA, PostGIS, …)?

I would suggest making tests with polygons of which you know the area (e.g. simple squares to begin with)

Moritz

–

James Duffy
PhD Researcher
Environment and Sustainability Institute
Penryn Campus
University of Exeter
Penryn
Cornwall
TR10 9FE

5

Just checked the data you send offline in QGIS 2.18.0.

With your sample data

SAGA (through processing) returns an area of 0.0000…

QGIS with OTF ON: 34,815 m²

QGIS with OTF OFF: 0,127 m²

Cheers,

Stefan

6

Le 19 mai 2017 10:19:06 GMT+02:00, "Blumentrath, Stefan" <Stefan.Blumentrath@nina.no> a écrit :

Just checked the data you send offline in QGIS 2.18.0.

With your sample data
SAGA (through processing) returns an area of 0.0000…
QGIS with OTF ON: 34,815 m²
QGIS with OTF OFF: 0,127 m²

In what projection is the original data ?

Moritz

Cheers,
Stefan

7

Output from ogrinfo:

INFO: Open of `C:\Sandbox\RGRASS\DATA\training_single.shp’

···

On 19 May 2017 at 09:38, Moritz Lennert <mlennert@club.worldonline.be> wrote:

Le 19 mai 2017 10:19:06 GMT+02:00, “Blumentrath, Stefan” <Stefan.Blumentrath@nina.no> a écrit :

Just checked the data you send offline in QGIS 2.18.0.

With your sample data
SAGA (through processing) returns an area of 0.0000…
QGIS with OTF ON: 34,815 m²
QGIS with OTF OFF: 0,127 m²

In what projection is the original data ?

Moritz

Cheers,
Stefan

–

James Duffy
PhD Researcher
Environment and Sustainability Institute
Penryn Campus
University of Exeter
Penryn
Cornwall
TR10 9FE

8

Le 19 mai 2017 10:41:29 GMT+02:00, James Duffy <james.philip.duffy@gmail.com> a écrit :

Output from ogrinfo:

INFO: Open of `C:\Sandbox\RGRASS\DATA\training_single.shp'

     using driver `ESRI Shapefile' successful.

Layer name: training_single

Metadata:

DBF_DATE_LAST_UPDATE=2017-05-17

Geometry: Polygon

Feature Count: 1

Extent: (73.213474, 0.225604) - (73.213477, 0.225608)

Layer SRS WKT:

GEOGCS["WGS 84",

   DATUM["WGS_1984",

       SPHEROID["WGS 84",6378137,298.257223563,

           AUTHORITY["EPSG","7030"]],

       AUTHORITY["EPSG","6326"]],

   PRIMEM["Greenwich",0,

       AUTHORITY["EPSG","8901"]],

   UNIT["degree",0.0174532925199433,

       AUTHORITY["EPSG","9122"]],

   AUTHORITY["EPSG","4326"]]

id: Integer64 (10.0)

Ok, so the data is in degrees and 0.127 is not in m2, but in degree2 which is not a really useful unit as its meaning on the ground depends on the latitude.

Moritz

On 19 May 2017 at 09:38, Moritz Lennert <mlennert@club.worldonline.be>
wrote:

Le 19 mai 2017 10:19:06 GMT+02:00, "Blumentrath, Stefan" <
Stefan.Blumentrath@nina.no> a écrit :
>Just checked the data you send offline in QGIS 2.18.0.
>
>With your sample data
>SAGA (through processing) returns an area of 0.0000…
>QGIS with OTF ON: 34,815 m²
>QGIS with OTF OFF: 0,127 m²

In what projection is the original data ?

Moritz

>
>Cheers,
>Stefan

9

On Thu, May 18, 2017 at 4:22 PM, James Duffy <james.philip.duffy@gmail.com> wrote:

Hello,

I’m trying to calculate the area of a polygon using v.to.db as follows (through R):

execGRASS(“v.to.db”, flags = c(“p”,“quiet”), parameters = list(map = “test”, units = “meters”, option = “area”, separator = “comma”),intern = TRUE)

I assume that the output of this command will be a meters squared calculation of the area inside the polygon ‘test’.

However, in QGIS, I manually calculated the area the they do not match. QGIS (which I believe to be true) shows 0.128m2, and the output from v.to.db shows 22.42226m2.

This is a large difference. Is it possible that several polygons share the same category? In this case the entry in the attribute table is the sum of all polygons with that category.

I am working in WGS84.

GRASS uses geodesic distance and area calculation for latlon, I believe QGIS as well.

Markus M

Can anyone shed some light on how v.to.db calculates area and why it might be going wrong please?

Thank you

James

–
James Duffy
PhD Researcher
Environment and Sustainability Institute
Penryn Campus
University of Exeter
Penryn
Cornwall
TR10 9FE

grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user

10

Le 19 mai 2017 11:12:30 GMT+02:00, Markus Metz <markus.metz.giswork@gmail.com> a écrit :

On Thu, May 18, 2017 at 4:22 PM, James Duffy
<james.philip.duffy@gmail.com>
wrote:

Hello,

I'm trying to calculate the area of a polygon using v.to.db as

follows
(through R):

execGRASS("v.to.db", flags = c("p","quiet"), parameters = list(map =

"test", units = "meters", option = "area", separator = "comma"),intern

TRUE)

I assume that the output of this command will be a meters squared

calculation of the area inside the polygon 'test'.

However, in QGIS, I manually calculated the area the they do not

match.
QGIS (which I believe to be true) shows 0.128m2, and the output from
v.to.db shows 22.42226m2.

This is a large difference. Is it possible that several polygons share
the
same category? In this case the entry in the attribute table is the sum
of
all polygons with that category.

I am working in WGS84.

GRASS uses geodesic distance and area calculation for latlon, I believe
QGIS as well.

QGIS only uses geodesic when OTF projection is activated and a datum is selected in the project properties (done by default with otf active I think). Otherwise it uses map units.

Moritz

Markus M

Can anyone shed some light on how v.to.db calculates area and why it

might be going wrong please?

Thank you

James

--
James Duffy
PhD Researcher
Environment and Sustainability Institute
Penryn Campus
University of Exeter
Penryn
Cornwall
TR10 9FE

_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user

11

Just to clear a few things up.

The shapefile contains only a single polygon (I created it as a simple test to work out this discrepancy).

OTF is turned off in QGIS, I am working in EPSG 4326.

The 0.127m2 comes from using the Measure tool (OTF off) in QGIS. And it is set to display m2.

Using the $area calculator in the field calculator in QGIS does indeed show something that is more likely degrees squared.

···

On 19 May 2017 at 10:20, Moritz Lennert <mlennert@club.worldonline.be> wrote:

Le 19 mai 2017 11:12:30 GMT+02:00, Markus Metz <markus.metz.giswork@gmail.com> a écrit :

On Thu, May 18, 2017 at 4:22 PM, James Duffy
<james.philip.duffy@gmail.com>
wrote:

Hello,

I’m trying to calculate the area of a polygon using v.to.db as
follows
(through R):

execGRASS(“v.to.db”, flags = c(“p”,“quiet”), parameters = list(map =
“test”, units = “meters”, option = “area”, separator = “comma”),intern
=
TRUE)

I assume that the output of this command will be a meters squared
calculation of the area inside the polygon ‘test’.

However, in QGIS, I manually calculated the area the they do not
match.
QGIS (which I believe to be true) shows 0.128m2, and the output from
v.to.db shows 22.42226m2.

This is a large difference. Is it possible that several polygons share
the
same category? In this case the entry in the attribute table is the sum
of
all polygons with that category.

I am working in WGS84.

GRASS uses geodesic distance and area calculation for latlon, I believe
QGIS as well.

QGIS only uses geodesic when OTF projection is activated and a datum is selected in the project properties (done by default with otf active I think). Otherwise it uses map units.

Moritz

Markus M

Can anyone shed some light on how v.to.db calculates area and why it
might be going wrong please?

Thank you

James

–
James Duffy
PhD Researcher
Environment and Sustainability Institute
Penryn Campus
University of Exeter
Penryn
Cornwall
TR10 9FE

grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user

–

James Duffy
PhD Researcher
Environment and Sustainability Institute
Penryn Campus
University of Exeter
Penryn
Cornwall
TR10 9FE

12

On 19/05/17 11:25, James Duffy wrote:

Just to clear a few things up.

The shapefile contains only a single polygon (I created it as a simple
test to work out this discrepancy).

OTF is turned off in QGIS, I am working in EPSG 4326.

The 0.127m2 comes from using the Measure tool (OTF off) in QGIS. And it
is set to display m2.

Using the $area calculator in the field calculator in QGIS does indeed
show something that is more likely degrees squared.

Below some results from quick tests.

Summary: personally I would trust GRASS GIS a bit more, but I'm biaised...

***************************
Test 1:

GRASS trunk in NC demo location (EPSG 3358)

g.region n=200100 s=200000 w=500000 e=500100 -ap
v.in.region out=nc_testpolygon
v.to.db -p nc_testpolygon op=area --q
1|10000

GRASS trunk in EPSG 4326 location
v.proj location=nc_spm_08 mapset=user1 input=nc_testpolygon
v.to.db -p nc_testpolygon op=area --q
1|10002.503642603
v.out.ogr nc_testpolygon out=nc_testpolygon.shp (i.e. shapefile is in degress, EPSG 4326)

QGIS (using the exported shapefile)

No OTF
measuring tool: 12.324,737 m²
$area in calculator: 12320.2346026

OTF to EPSG 4326, with datum WGS84
measuring tool: 10.072,355 m²
$area in calculator: 10046.3679510

********************************
Test 2:

GRASS trunk, in ESPG 4326 location

g.region n=74 s=73 w=0 e=1 res=1 -p
v.in.region region_1degree_73_74N
v.to.db -p region_1degree_73_74N op=area --q
1|3539348689.40454

v.out.ogr region_1degree_73_74N out=region_1degree_73_74N.shp

GRASS trunk in EPSG 3575 location
v.proj location=LL_WGS84 mapset=mlennert input=region_1degree_73_74N
v.to.db -p region_1degree_73_74N op=area --q
1|3539337483.97153

QGIS (using the exported shapefile)
no OTF projection

measuring tool 1,002 m² (sic)
export geometry column: 1
$area in calculator 12392029030.5000000

with OTF projection, datum for measurements WGS84

measuring tool 3.503.322.774,337 m²
export geometry column: 1
$area in calculator: 3499702798.3400998

reproject to EPSG 3575
with OTF, datum WGS84

export geometry column 3539169026.18628
export geoemtry column (column named otf) 3499702798.00000
$area in calculator: 3499702798.6579294

13

Hello,

Thank you for your investigations.

Is there anything else I can do about this? I trust my QGIS measurement because I know the size of the object I have drawn around. However I would really like to use GRASS to automate my workflow…

Thanks

James

···

On 19 May 2017 at 13:41, Moritz Lennert <mlennert@club.worldonline.be> wrote:

On 19/05/17 11:25, James Duffy wrote:

Just to clear a few things up.

The shapefile contains only a single polygon (I created it as a simple
test to work out this discrepancy).

OTF is turned off in QGIS, I am working in EPSG 4326.

The 0.127m2 comes from using the Measure tool (OTF off) in QGIS. And it
is set to display m2.

Using the $area calculator in the field calculator in QGIS does indeed
show something that is more likely degrees squared.

Below some results from quick tests.

Summary: personally I would trust GRASS GIS a bit more, but I’m biaised…

Test 1:

GRASS trunk in NC demo location (EPSG 3358)

g.region n=200100 s=200000 w=500000 e=500100 -ap
v.in.region out=nc_testpolygon
v.to.db -p nc_testpolygon op=area --q
1|10000

GRASS trunk in EPSG 4326 location
v.proj location=nc_spm_08 mapset=user1 input=nc_testpolygon
v.to.db -p nc_testpolygon op=area --q
1|10002.503642603
v.out.ogr nc_testpolygon out=nc_testpolygon.shp (i.e. shapefile is in degress, EPSG 4326)

QGIS (using the exported shapefile)

No OTF
measuring tool: 12.324,737 m²
$area in calculator: 12320.2346026

OTF to EPSG 4326, with datum WGS84
measuring tool: 10.072,355 m²
$area in calculator: 10046.3679510

Test 2:

GRASS trunk, in ESPG 4326 location

g.region n=74 s=73 w=0 e=1 res=1 -p
v.in.region region_1degree_73_74N
v.to.db -p region_1degree_73_74N op=area --q
1|3539348689.40454

v.out.ogr region_1degree_73_74N out=region_1degree_73_74N.shp

GRASS trunk in EPSG 3575 location
v.proj location=LL_WGS84 mapset=mlennert input=region_1degree_73_74N
v.to.db -p region_1degree_73_74N op=area --q
1|3539337483.97153

QGIS (using the exported shapefile)
no OTF projection

measuring tool 1,002 m² (sic)
export geometry column: 1
$area in calculator 12392029030.5000000

with OTF projection, datum for measurements WGS84

measuring tool 3.503.322.774,337 m²
export geometry column: 1
$area in calculator: 3499702798.3400998

reproject to EPSG 3575
with OTF, datum WGS84

export geometry column 3539169026.18628
export geoemtry column (column named otf) 3499702798.00000
$area in calculator: 3499702798.6579294

–

James Duffy
PhD Researcher
Environment and Sustainability Institute
Penryn Campus
University of Exeter
Penryn
Cornwall
TR10 9FE

14

On 24/05/17 11:06, James Duffy wrote:

Hello,

Thank you for your investigations.

Is there anything else I can do about this? I trust my QGIS measurement
because I know the size of the object I have drawn around.

But as you say below:

        OTF is turned off in QGIS, I am working in EPSG 4326.

        The 0.127m2 comes from using the Measure tool (OTF off) in QGIS.
        And it
        is set to display m2.

As I showed with my test, the value of 0.127m2 is probably wrong. In the case of the 10000m2 rectangle, I get an error of over 2000m2 with OTF off on an EPSG 4326 ! Only with OTF on, I get a result closer to the true size.

Would you be willing to share your test polygon ?

Moritz

However I
would really like to use GRASS to automate my workflow...

Thanks

James

On 19 May 2017 at 13:41, Moritz Lennert <mlennert@club.worldonline.be
<mailto:mlennert@club.worldonline.be>> wrote:

    On 19/05/17 11:25, James Duffy wrote:

        Just to clear a few things up.

        The shapefile contains only a single polygon (I created it as a
        simple
        test to work out this discrepancy).

        OTF is turned off in QGIS, I am working in EPSG 4326.

        The 0.127m2 comes from using the Measure tool (OTF off) in QGIS.
        And it
        is set to display m2.

        Using the $area calculator in the field calculator in QGIS does
        indeed
        show something that is more likely degrees squared.

    Below some results from quick tests.

    Summary: personally I would trust GRASS GIS a bit more, but I'm
    biaised...

    ***************************
    Test 1:

    GRASS trunk in NC demo location (EPSG 3358)

    g.region n=200100 s=200000 w=500000 e=500100 -ap
    v.in.region out=nc_testpolygon
    v.to.db -p nc_testpolygon op=area --q
    1|10000

    GRASS trunk in EPSG 4326 location
    v.proj location=nc_spm_08 mapset=user1 input=nc_testpolygon
    v.to.db -p nc_testpolygon op=area --q
    1|10002.503642603
    v.out.ogr nc_testpolygon out=nc_testpolygon.shp (i.e. shapefile is
    in degress, EPSG 4326)

    QGIS (using the exported shapefile)

    No OTF
    measuring tool: 12.324,737 m²
    $area in calculator: 12320.2346026

    OTF to EPSG 4326, with datum WGS84
    measuring tool: 10.072,355 m²
    $area in calculator: 10046.3679510

    ********************************
    Test 2:

    GRASS trunk, in ESPG 4326 location

    g.region n=74 s=73 w=0 e=1 res=1 -p
    v.in.region region_1degree_73_74N
    v.to.db -p region_1degree_73_74N op=area --q
    1|3539348689.40454

    v.out.ogr region_1degree_73_74N out=region_1degree_73_74N.shp

    GRASS trunk in EPSG 3575 location
    v.proj location=LL_WGS84 mapset=mlennert input=region_1degree_73_74N
    v.to.db -p region_1degree_73_74N op=area --q
    1|3539337483.97153

    QGIS (using the exported shapefile)
    no OTF projection

    measuring tool 1,002 m² (sic)
    export geometry column: 1
    $area in calculator 12392029030 <tel:12392029030>.5000000

    with OTF projection, datum for measurements WGS84

    measuring tool 3.503.322.774,337 m²
    export geometry column: 1
    $area in calculator: 3499702798.3400998

    reproject to EPSG 3575
    with OTF, datum WGS84

    export geometry column 3539169026.18628
    export geoemtry column (column named otf) 3499702798.00000
    $area in calculator: 3499702798.6579294

--
*James Duffy*
PhD Researcher
Environment and Sustainability Institute
Penryn Campus
University of Exeter
Penryn
Cornwall
TR10 9FE