area not calculated correctly with OTF on #21270

qgib · 2015-08-11T12:09:46Z

Author Name: Andrew McAninch (@spaceof7)
Original Redmine Issue: 13209
Affected QGIS version: master
Redmine category:field_calculator
Assignee: Nyall Dawson

I had problems calculating areas in 2.8.3 and found several related reported bugs that were listed as closed and fixed in Master. I tested in Master but I'm getting errors.

I have a couple datasets that were created by buffering stream lines with CRS EPSG:2285. With OTF turned off the areas are the same as those calculated in spatialite and arcgis. If OTF is on(same results with either the same CRS as the layer or different CRS) and the ellipsoid is left as the default GRS 1980 the calculated areas are way off, 22,365 vs 241,226. I get the same results whether the units are set to meters or feet. finally, if i set the ellipsoid to planimetric then the areas are calulated as 0.

buffer_dnr.zip (Andrew McAninch)
knox.zip (Randal Hale)
area_test.zip (Andrew McAninch)
area_test_3857.cpg.zip (Giovanni Manghi)
10000sqft-test.zip (Steven Mizuno)
QGIS-_13209-testing-results.xlsx (Steven Mizuno)

Related issue(s): #18255 (relates), #20259 (relates), #22092 (relates)
Redmine related issue(s): 9690, 12057, 14082

qgib · 2015-10-09T13:53:00Z

Author Name: John Tull (John Tull)

I can confirm that this bug exists in trunk as well. With OTF on, area calculations are many orders of magnitude off for data in a UTM meters coordinate system. Areas are accurately calculated when OTF is turned off.

qgib · 2015-10-10T05:20:48Z

Author Name: Randal Hale (@rjhale1971)

There are several of these bug reports floating around #18255 -> this is one that is I think still open.

qgib · 2015-10-10T05:21:39Z

Author Name: Randal Hale (@rjhale1971)

fixed_version_id was configured as Version 2.12

qgib · 2015-10-10T05:22:51Z

Author Name: Randal Hale (@rjhale1971)

I followed #18255 and it's still incorrect - although in Master (2.11) area gets set to 0. Which is slightly better than something that looks slightly right.

qgib · 2015-10-10T05:32:13Z

Author Name: Randal Hale (@rjhale1971)

Seems to also be related to this one that is closed #20259

qgib · 2015-10-14T22:49:30Z

Author Name: Nyall Dawson (@nyalldawson)

Please test with current master, this may have been fixed along with #21643

status_id was changed from Open to Feedback

qgib · 2015-10-15T09:31:42Z

Author Name: Andrew McAninch (@spaceof7)

I just tested with the latest OSGEO4W nightly(2.11.0-88) and I get the same errors.

qgib · 2015-10-15T10:21:14Z

Author Name: Nyall Dawson (@nyalldawson)

You'll need the next build, that one is too early

qgib · 2015-10-16T08:58:14Z

Author Name: Andrew McAninch (@spaceof7)

I updated to 2.11.0-89 and the bug is not completely fixed. Now if OTF is on and the ellipsoid is set the "None/Planimetric" the areas are calculated correctly. However, if OTF is on and I leave the ellipsoid set to the default then the values are wrong.

qgib · 2015-10-16T16:12:12Z

Author Name: Nyall Dawson (@nyalldawson)

assigned_to_id was configured as Nyall Dawson

qgib · 2015-10-16T19:31:10Z

Author Name: Nyall Dawson (@nyalldawson)

Hmm... I can't reproduce this, using the attached dataset and described steps. Here's what I see:

For the feature which sqft_arc = '142434'

OTF OFF: area = 1.323 ha
OTF ON, EPSG 2285, default ellipsoid (GRS 1980): area = 1.322 ha
OTF ON, EPSG 2285, ellipsoid = none/planimetric: area = 1.323 ha
OTF ON, EPSG 4326, ellipsoid = WGS84: area = 1.322 ha

All looks OK to me. This is using the identify tool and looking at the derived area attribute. What method are you using to determine the area?

qgib · 2015-10-17T03:17:16Z

Author Name: Andrew McAninch (@spaceof7)

I should have specified that I am using the field calculator. I tested the identify tool and i do get the correct areas using that, its just when I calc a new field that I get wrong values.

qgib · 2015-10-17T10:46:30Z

Author Name: Nyall Dawson (@nyalldawson)

Is this using a virtual field? Could it be a duplicate of #20739?

qgib · 2015-10-17T13:23:02Z

Author Name: Stefan Blumentrath (@ninsbl)

Hi Andrew, QGIS computes - different from almost all other GIS - by default an ellipsoidal area measurement, and not one a kartesian one (at least as long as you do not set the measurement ellipsoid manually to "None/Planimetric" or turn off OTF-projection). That means e.g. if you digitize a square with 1000 x 1000 m boundary length, the returned area for it will not be 1 km2 but a slightly different number (e.g. 0.99986789 km2), depending on where you square is located relatively to the central median. You can find a quite long discussion about it here: #20259
So, unless your results are not completely off, the difference you (and also Nyall - even if it is just 0.001ha in this case) are reporting may be intended behavior...

qgib · 2015-10-17T15:35:51Z

Author Name: Andrew McAninch (@spaceof7)

Nyall - It happens when I am adding a new field, not a virtual one, well, I haven't tested it with a virtual field, so it may behave differently with a virtual field. I can test that as well.

Stefan - Thanks, I am aware of this, but the errors I get are the same as I mentioned in the original bug report, which is to say they are off by at least an order of magnitude.

qgib · 2015-10-17T15:54:32Z

Author Name: Andrew McAninch (@spaceof7)

I tested this with a virtual field noticed something unexpected. I am looking at a feature where the correct area should be 3,988,722. With OTF on, when I calculate a new virtual field the output preview in the field calculator window is the same incorrect value I get when adding a new field, 369,889. However, the value that actually gets filled in the virtual field is 3988722, the correct value.

qgib · 2015-10-17T16:32:53Z

Author Name: Randal Hale (@rjhale1971)

I'm reading through All of this and I have a practical example:

Attached is a shapefile in EPSG:2240. There is an area field added. Calculate the area with OTF on and OTF off using the field calculator. The results are wildly varying. Converting the results to Acres gives me 4.7(ish) with OTF on and 50.5 with OTF off. The area is 50.5 acres (or 2200346 square feet). This has been an on going problem on my end starting at 2.4 or 2.6. It's not ellipsoidal. I haven't tested this on windows - I'm on Ubuntu 14.04.3 2.10 QGIS.

It seems like it's returning almost sq meters although the projection is in US Feet. Which almost maybe the ellipsoidal issue mentioned.

Anyway - I always hope this is just me - but if it is - possibly more of an indication of insanity. I don' tthink it is - but I will never rule out insanity.

9205 was configured as knox.zip

knox.zip (Randal Hale) - Shapefile in EPSG2240

qgib · 2015-10-17T22:03:45Z

Author Name: Mathieu Pellerin - nIRV (Mathieu Pellerin - nIRV)

Playing around with Randal's shapefile, here are a couple of observations:

the area shown (via identify tool) when OTF is off is 20.442 ha, and when OTF is on it is 20.538 ha which amounts to a 0.096 ha difference (0.5% more or less) and I have no idea if that is an acceptable variance or not
when OTF is off, the field calculator returns the proper value for $area, that is 2200346.6519 feet
when OTF is on, the field calculator is definitively broken (whether you create a new field, or update an existing field is irrelevant), and the $area value returned through it is 205380.4847
interestingly, if you label the polygon using the "$area" expression, that value is correct irrespective of whether OTF is on or off

qgib · 2015-10-18T02:00:23Z

Author Name: Nyall Dawson (@nyalldawson)

Hmmm... I can't work out what the expected behaviour is here, and have run out of time for 2.12 fixes.

assigned_to_id removed Nyall Dawson

qgib · 2015-10-18T03:52:35Z

Author Name: Randal Hale (@rjhale1971)

This is one of those things that keeps getting lumped into the Ellipsoid issue. Heh. Which is quite maddening once you see the long winded discussion and then this Field Calculator thing just disappears into the fog. No worries - I know 2.12 was getting close so maybe this will make it during 2.13 for 2.14. As long as someone at the top knows it's happening I feel better.

qgib · 2015-10-18T17:00:46Z

Author Name: Mathieu Pellerin - nIRV (Mathieu Pellerin - nIRV)

Randal,

What I see here is that when OTF is on, the $area value is returned in square meters (i.e. 205380.4847 = 20.538ha) which is a near-perfect match to the square foot value return when OTF is off (that is 2200346.6519 = 20.442ha).

Can you confirm this is the case with you too? If not, can you write down a precise step to reproduce any value that doesn't match the above?

qgib · 2015-10-18T17:23:12Z

Author Name: Randal Hale (@rjhale1971)

So when I go in (and check my math please - I've been working all day)

I've turned File -> Project Properties -> General -> Measure Tool -> Ellipsoid to Non/Planimetric (just to cover myself)

A. If I calculate the area with OTF off I get 2200346.65259 with EPSG:2240
B. OTF On I get 205380.48474 with EPSG:2240

If my math is correct to convert Sq Feet to Sq Meters I would multiply the value in A * 0.092903 (I blame google if this is wrong - ha) and I get 204418.805065569 (which is close)

If I project my shapefile to EPSG:26916 I get 204340.82666 with OTF off and 205381.07575 with OTF on.....

So I would say yes - you are correct it's square meters. Apparently whatever is happening isn't happening when I'm using a projection in meters (or it's happening very very little).

If that makes sense.

qgib · 2015-11-07T10:26:02Z

Author Name: Giovanni Manghi (@gioman)

fixed_version_id removed Version 2.12

qgib · 2015-11-10T02:24:41Z

Author Name: Giovanni Manghi (@gioman)

It seems to me that in QGIS master things are ok. There is just still an issue with virtual fields

#20739

category_id was configured as Field calculator

qgib · 2015-11-12T12:25:29Z

Author Name: Andrew McAninch (@spaceof7)

Giovanni,
I am still getting the same behaviour when calculating a non-virtual field with 2.13-16.

Mathieu,
That seems to be the case(with OTF on $area is calculated in square meters) for MOST of the features in the dataset I am testing but there are a few(11 out of 638) where the areas are WAY off. For all features the areas are calculated correctly when OTF is off and also with OTF on and the ellipsoid is set to non/planimetric. but with OTF on and the default ellipsoid I get bad values, for instance:

correct area(m2) calculated area(m2)
9,213 12,829,041
13,437 3,686,013
140,552 16,559,293
855,174 559,466

The steps I took were:
field calculator button
new field, integer = $area

Im attaching the shapefile I am using to test this with several test fields calculated so you can see the values I am getting.

9313 was configured as area_test.zip

area_test.zip (Andrew McAninch)

qgib · 2015-11-13T09:23:09Z

Author Name: Giovanni Manghi (@gioman)

Andrew McAninch wrote:

Giovanni,
I am still getting the same behaviour when calculating a non-virtual field with 2.13-16.

JC, I see... To be more clear I attached here your dataset but in 3857, to units are meters (because we know there is also the issue of QGIS computing always in meters even if the CRS is in feet). I left 3 columns, area by postgis, area by qgis no OTF, area by qgis with OTF. I have tested also other desktop gis packages and the results are correct.

priority_id was changed from Normal to Severe/Regression
9315 was configured as area_test_3857.cpg.zip
status_id was changed from Feedback to Open

area_test_3857.cpg.zip (Giovanni Manghi)

qgib · 2015-11-27T10:55:48Z

Author Name: Steven Mizuno (Steven Mizuno)

The problems noted in this ticket (and some others related to it) interested me as I had been considering how to go about testing QGIS. This seemed to be an excellent opportunity to refine some procedures and hopefully define the problems coherently.

I have concentrated on the data in area_test for my testing. One thing I noticed in this data is that the various calculated areas are integer values. This made it hard to spot small differences, so I have used floating point storage for my testing.

And I got curious about what happens if the geometries are changed -- moving or altering the shape.

I transformed area_test from EPSG:2285 to EPSG:32148 to have the data in the same basic projection, but with meters for units. This is easily done with PostGIS or SpatiaLite using DB Manager to import/export.
I considered that transforming to EPSG:3857, as some have done, throws in complications for finding the root cause of the problems.

The 10000sqft-test-2285 file was created to have a very specific known area for some tests.

Note that in the procedures below, when I indicate "OTF on" this means only checking "Enable 'on the fly' CRS transformation" in Project Properties and accepting the ellipsoid that the CRS uses, in this case, GRS 1980. I did not test the effect of using different ellipsoids.

Not mentioned before, but there is another function that calculates area: [[area()]] -- more specifically, [[area($geometry)]] would appear to be equal to [[$area]] (which it isn't.)
I'm not sure whether area() is strictly plane or can handle ellipsoidal calculations.

Also not mentioned is any possible effect of "Preferred units" in Settings|Options|Map Tools, so I have tested this. The only effect I saw was in Identify Results.

I haven't included much on what Identify Results, Derived Area shows, but generally it does show the same area as calculated by $area, but the value is rounded so it is hard to be sure.

Based on testing I conducted:

A. when OTF is on, the $area calculation has a problem with geometries. Specifically, the area calculation based on an ellipsoid is defective. A few geometries in the provided area_test file are very obviously incorrect, but all aren't really close to what they should be.

B. two different area calculations are used: one for [[$area]]; a different one for [[area($geometry)]] -- the results are very slightly different for planar area calculation. This is not obvious on geometries that have whole numbers for vertex x, y values and have regular shapes like rectangles.

This leads to confusion because of the very slightly different results.

C. from running some PostGIS queries, I believe that all of the QGIS area computations on area_test using [[$area]] with OTF on are likely incorrect by a significant amount. The area computed based on ellipsoid should be within less than 0.1% of the value in sqm_arc field, I should think, as the CRS is a State Plane Coordinates projection.

D. the Field calculator has confusing behavior on actual field vs. virtual field

Field Calculator: the $area value placed in a virtual field uses the planar area value instead of the area based on ellipsoid; this is unexpected
Field Calculator: the $area value placed in an actual field is the area based on ellipsoid in meters

E. when OTF is on, in Expression Builder, for a virtual field, the preview value for $area varies with whether the builder is used to create the field or is editing an existing field. This is very confusing.

data showing the above points:
[using QGIS dev 5b59609]

A.

area_test feature ID=301/sqft_arc=144636: $area = 144636.18460 sq ft (OTF off) and 5626033.17938 (OTF on)

These two were very interesting and shouldn't happen (using feature ID= 301 and 303):

changing the shape of specific objects -- does not appreciably change the area calculated IF the area was way off.
moving specific objects -- changes the area calculated IF the area calc was way off
tested the calculated area: OTF to non-OTF and found that all features were different by about 0.1% or more; at first, this seemed reasonable (except for a few large differences), but on reconsideration I would expect a smaller difference as the CRS is a State Plane Coordinates projection (see PostGIS test below)
I tested both area_test and area_test transformed to EPSG:32148 Washington N in meters, with similar results for both

------------
testing the percent change of QGIS calculated area: OTF area to non-OTF area
Using Field Calculator create the following fields/expressions on the sample data (use type double, 19 / 10 ):

actual field: area_noOTF	$area *((1200.0/3937.0)^2)	(NOTE: turn off OTF first; convert to sq m -- don't need to convert if units is m)
actual field: area_OTF		$area				(NOTE: turn on OTF first)
virtual or actual field: pct_chg	("area_OTF"/("area_noOTF")*100)-100 

------------

B. example: for the area_test feature with sqft_arc=144636: $area=144636.18460083, area($geometry)=144636.18456012; from PostGIS: st_area()=144636.18456012, which is the same as area($geometry)

C. after importing area_test to PostGIS I ran ran some area computations; following is the query I used to compare the area determined from geography type (WGS84 ellipsoid) to the area of the original plane projection geometry type. Note that geography type is limited to the WGS84 ellipsoid, but has no significant difference compared to GRS80.

------------
-- PostGIS 2.1  to get area in various CRS (SRS in PostGIS parlance); compare 2 particular results
select id, area_2285_m, area_32148, area_26910, area_4326_geog,
	( area_4326_geog / area_2285_m * 100 ) - 100 as pct_chg	-- compare geography area (sq m) to cartesian area (sq ft conv. to sq m)
from 
(
	select id, 
		st_area(geom) * ((1200.0/3937.0)^2) as area_2285_m,	-- Washington N, ft; planar ft->m
		st_area(st_transform(geom,32148)) as area_32148,	-- Washington N, m; planar
		st_area(st_transform(geom,26910)) as area_26910,	-- UTM zone 10; planar
		st_area(st_transform(geom,4326)::geography) as area_4326_geog	-- ellipsoid, m
	from test_data.area_test
) x
order by id;
------------

The pct_chg is not more than approx. +/- 0.025% (about 2.5 parts in 10000)

Note that State Plane Coordinates (2285 & 32148) are intended to be accurate to about 1 part in 10000 for position, as I understand it. Area accuracy would be somewhat less.

As I don't know what the accuracy of the PostGIS (Proj4 and GEOS internally) functions are, I believe the PostGIS computed areas to be reasonable.

Projections like UTM cover large areas so they tend to be less accurate.

D. 1. & D. 2. see Excel worksheet: QGIS-#21270-testing-results.xlsx
this shows changing three variables and noting the results for the 10000sqft test files

E. observe: on virtual field creation, for $area the preview value is the calculated ellipsoid value in meters; on editing the expression the preview is the planar value in layer units.

I have chosen not to include the results as anyone should be able to reproduce the results from the procedures provided. I have noted the specific QGIS development version I used as a different one may produce different results, depending on what has been changed.

I am including an Excel worksheet showing the results for point D. as it is rather tedious to construct.

And I am also including the two files used for point D. containing a 10000ft x 1ft test object in CRS EPSG:2285 and EPSG:32148
10000sqft-test.zip containing: 10000sqft-test-2285.shp and 10000sqft-test-32148.shp

I believe the tests I have put forth show that the following things need to be fixed:

the ellipsoidal area calculation used by [[$area]] and Identify (most important!)
for plane area calculation [[$area]] should be equal to [[area($geometry)]]
the logic of Field Calculator as to how [[$area]] is calculated and shown as well as stored in a field (actual vs. virtual) should be consistent

9359 was configured as QGIS-_13209-testing-results.xlsx
9358 was configured as 10000sqft-test.zip

QGIS-_13209-testing-results.xlsx (Steven Mizuno) - results for point D (Excel worksheet)
10000sqft-test.zip (Steven Mizuno) - 10000ft x 1ft test feature in 2 CRS

qgib · 2016-01-06T04:07:03Z

Author Name: Giovanni Manghi (@gioman)

area not calculated correctly with OTF on #21270

area not calculated correctly with OTF on #21270

Comments

qgib commented Aug 11, 2015 • edited Loading

qgib commented Oct 9, 2015

qgib commented Oct 10, 2015 • edited Loading

qgib commented Oct 10, 2015

qgib commented Oct 10, 2015 • edited Loading

qgib commented Oct 10, 2015 • edited Loading

qgib commented Oct 14, 2015 • edited Loading

qgib commented Oct 15, 2015

qgib commented Oct 15, 2015

qgib commented Oct 16, 2015

qgib commented Oct 16, 2015

qgib commented Oct 16, 2015

qgib commented Oct 17, 2015

qgib commented Oct 17, 2015 • edited Loading

qgib commented Oct 17, 2015 • edited Loading

qgib commented Oct 17, 2015

qgib commented Oct 17, 2015

qgib commented Oct 17, 2015

qgib commented Oct 17, 2015

qgib commented Oct 18, 2015

qgib commented Oct 18, 2015

qgib commented Oct 18, 2015

qgib commented Oct 18, 2015

qgib commented Nov 7, 2015

qgib commented Nov 10, 2015 • edited Loading

qgib commented Nov 12, 2015

qgib commented Nov 13, 2015

qgib commented Nov 27, 2015 • edited Loading

qgib commented Jan 6, 2016 • edited Loading

qgib commented Feb 9, 2016

qgib commented Feb 14, 2016 • edited Loading

qgib commented Feb 16, 2016

qgib commented Feb 21, 2016

qgib commented Feb 21, 2016

qgib commented Feb 21, 2016

qgib commented Feb 21, 2016

qgib commented Feb 21, 2016

qgib commented Feb 25, 2016

qgib commented Aug 11, 2015 •

edited

Loading

qgib commented Oct 10, 2015 •

edited

Loading

qgib commented Oct 10, 2015 •

edited

Loading

qgib commented Oct 10, 2015 •

edited

Loading

qgib commented Oct 14, 2015 •

edited

Loading

qgib commented Oct 17, 2015 •

edited

Loading

qgib commented Oct 17, 2015 •

edited

Loading

qgib commented Nov 10, 2015 •

edited

Loading

qgib commented Nov 27, 2015 •

edited

Loading

qgib commented Jan 6, 2016 •

edited

Loading

qgib commented Feb 14, 2016 •

edited

Loading