Bug report #13209

area not calculated correctly with OTF on

Added by Andrew McAninch almost 9 years ago. Updated over 8 years ago.

Status:Closed
Priority:Severe/Regression
Assignee:Nyall Dawson
Category:Field calculator
Affected QGIS version:master Regression?:No
Operating System: Easy fix?:No
Pull Request or Patch supplied:No Resolution:fixed/implemented
Crashes QGIS or corrupts data:No Copied to github as #:21270

Description

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 (1.95 MB) Andrew McAninch, 2015-08-11 02:09 PM

knox.zip - Shapefile in EPSG2240 (4.42 KB) Randal Hale, 2015-10-17 06:32 PM

area_test.zip (1.95 MB) Andrew McAninch, 2015-11-12 01:25 PM

area_test_3857.cpg.zip (1.8 MB) Giovanni Manghi, 2015-11-13 10:23 AM

10000sqft-test.zip - 10000ft x 1ft test feature in 2 CRS (3.47 KB) Steven Mizuno, 2015-11-27 11:55 AM

QGIS-_13209-testing-results.xlsx - results for point D (Excel worksheet) (11.3 KB) Steven Mizuno, 2015-11-27 11:55 AM


Related issues

Related to QGIS Application - Bug report #12057: Computed area is wrong when reprojection is active Closed 2015-01-26
Related to QGIS Application - Bug report #9690: Incorrect areas using field calculator. Closed 2014-03-02
Related to QGIS Application - Bug report #14082: area calculation in Field Calculator is depending on Outp... Closed 2016-01-10

Associated revisions

Revision 229ef291
Added by Nyall Dawson over 8 years ago

If no geomCalculator is set for QgsExpression, perform cartesian
calculations for $length, $area, $perimeter (refs #13209)

ie, the default is now to use planimeteric calculations unless
a geomCalculator has been explicitly set

Revision 9f772bce
Added by Nyall Dawson over 8 years ago

Add API call for setting preferred distance unit for QgsExpression

(refs #13209, #12939, #4857)

Revision 6453907c
Added by Nyall Dawson over 8 years ago

Fix project unit confusion (part 1)

- Remove existing (confusing and totally useless) "canvas units"
setting. This setting was supposed to be used for a confusing number
of totally different uses, eg changing the coordinate display
format and changing the default measurement units, but it was
so broken it did none of these things except change the
coordinate display and it only did that if OTF was off.

- Add a new "coordinate display" section in project settings,
which allows choice of format for coordinate display via a
combo box (this will also make it easy to add additional custom
formats in future), and make this setting work regardless of
whether OTF is on or off. This setting applies to both the
coordinate display in the status bar and coordinates shown
via the identify tool

(refs #13209, fixes #9730)

Revision 17a29f90
Added by Nyall Dawson over 8 years ago

Add area units to QgsUnitTypes (refs #13209)

Revision ddbdcf8a
Added by Nyall Dawson over 8 years ago

Fix project unit confusion (pt 2): add project distance unit setting

Adds a new option in project properties to set the units used for
distance measurements. This setting defaults to the units set in
QGIS options, but can then be overridden for specific projects.

The setting is respected for length and perimeter calculations in:
- Attribute table field update bar
- Field calculator calculations
- Identify tool derived length and perimeter values
- Default unit shown in measure dialog

Also adds unit tests to ensure that length and perimeter calculated
by attribute table update bar, field calculator and identify tool
are consistent wrt ellipsoidal calculations and distance units.

(refs #13209, #12939, #2402, #4857, #4252)

Revision dfdcec89
Added by Nyall Dawson over 8 years ago

Fix project unit confusion (pt 3): add area unit settings with a
ton of available area units (eg m2, km2, mi2, ft2, yd2, ha, ac,
etc)

Adds a new option in both the QGIS setting and project properties to
set the units used for area measurements. Just like the distance
setting, this defaults to the units set in QGIS options, but can
then be overridden for specific projects.

The setting is respected for area calculations in:
- Attribute table field update bar
- Field calculator calculations
- Identify tool derived length and perimeter values

Also adds unit tests to ensure that area calculated by attribute table
update bar, field calculator and identify tool are consistent wrt
ellipsoidal calculations and area units.

TODO: make measure tool respect area setting

(refs #13209, #4252 and fixes #12939, #2402, #4857)

Revision 479d90a5
Added by Nyall Dawson over 8 years ago

Make measure tool respect project area unit setting

Also add unit tests for measure tool length and area measurement
to ensure they return the same results as field calculator
and identify tool

Refs #13209, #4252

History

#1 Updated by John Tull over 8 years ago

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.

#2 Updated by Randal Hale over 8 years ago

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

#3 Updated by Randal Hale over 8 years ago

  • Target version set to Version 2.12

#4 Updated by Randal Hale over 8 years ago

I followed #9690 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.

#5 Updated by Randal Hale over 8 years ago

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

#6 Updated by Nyall Dawson over 8 years ago

  • Status changed from Open to Feedback

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

#7 Updated by Andrew McAninch over 8 years ago

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

#8 Updated by Nyall Dawson over 8 years ago

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

#9 Updated by Andrew McAninch over 8 years ago

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.

#10 Updated by Nyall Dawson over 8 years ago

  • Assignee set to Nyall Dawson

#11 Updated by Nyall Dawson over 8 years ago

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?

#12 Updated by Andrew McAninch over 8 years ago

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.

#13 Updated by Nyall Dawson over 8 years ago

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

#14 Updated by Stefan Blumentrath over 8 years ago

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: #12057
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...

#15 Updated by Andrew McAninch over 8 years ago

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.

#16 Updated by Andrew McAninch over 8 years ago

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.

#17 Updated by Randal Hale over 8 years ago

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.

#18 Updated by Mathieu Pellerin - nIRV over 8 years ago

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

#19 Updated by Nyall Dawson over 8 years ago

  • Assignee deleted (Nyall Dawson)

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

#20 Updated by Randal Hale over 8 years ago

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.

#21 Updated by Mathieu Pellerin - nIRV over 8 years ago

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?

#22 Updated by Randal Hale over 8 years ago

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.

#23 Updated by Giovanni Manghi over 8 years ago

  • Target version deleted (Version 2.12)

#24 Updated by Giovanni Manghi over 8 years ago

  • Category set to Field calculator

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

#12622

#25 Updated by Andrew McAninch over 8 years ago

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.

#26 Updated by Giovanni Manghi over 8 years ago

  • Priority changed from Normal to Severe/Regression
  • File area_test_3857.cpg.zip added
  • Status changed from Feedback to Open

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.

#27 Updated by Steven Mizuno over 8 years ago

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
1. 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
2. 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-#13209-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

#28 Updated by Giovanni Manghi over 8 years ago

see also #14057

#29 Updated by Nyall Dawson over 8 years ago

  • Assignee set to Nyall Dawson
  • Status changed from Open to In Progress
  • % Done changed from 0 to 50

#30 Updated by Nyall Dawson over 8 years ago

  • Status changed from In Progress to Feedback
  • % Done changed from 50 to 100

With 479d90a5ceec1c1b658aa101702ecf855e534987 this should now be resolved, with the one exception of virtual fields (#12622).

I'd appreciate extensive testing and feedback so that we can get this closed before 2.14 is released.

#31 Updated by Andrew McAninch over 8 years ago

Excellent! I briefly tested this in master with the data I had been using and the areas are calculating correctly. I will try to do some more thorough testing of measurements this week.

#32 Updated by Randal Hale over 8 years ago

So I just loaded master on ubuntu (pulling from deb http://qgis.org/debian-nightly wily main) and it still appears to not be working correctly with OTF checked. I checked a polygon that was 834 acres and did a $area/43560 and it appears to still be defaulting to meters. My projection is EPSG:2274. I'm still checking it but I wanted to leave a comment.

#33 Updated by Nyall Dawson over 8 years ago

Randal - did you check the new setting under project properties for area units? What's that set to?

#34 Updated by Randal Hale over 8 years ago

Holy Crap - It was set to square meters. I just switched it to square feet and it's working. You are my hero. I didn't know that was an option.

#35 Updated by Nyall Dawson over 8 years ago

It's a new option, introduced to help fix this issue. You may also want to check the "default units" under qgis options - > map tools. That's where you set the default values for the units for new projects.

#36 Updated by Randal Hale over 8 years ago

I'm so happy I could jump up and down.

#37 Updated by Nyall Dawson over 8 years ago

  • Status changed from Feedback to Closed
  • Resolution set to fixed/implemented

Also available in: Atom PDF