4

I have a polygon shapefile in geographic coordinate system, I would like to calculate area in square meters in QGIS using Python.

Is there any formula or equation to do this?

5
  • 1
    Are you looking for a Python script or just the formula to calculate the area for each polygon? There's a $area formula used in the Field Calculator but I am unsure how to incorporate that into Python.
    – Joseph
    Commented Nov 12, 2014 at 14:05
  • 1
    Questions about measuring lengths/areas with data in geographic coordinates are very common nowadays. What you need to do is to project your data into a suitable projected coordinate system that is using meters as units. What is the suitable CRS is something that has been dealt with in gis.stackexchange.com/questions/121489/…
    – user30184
    Commented Nov 12, 2014 at 14:45
  • I required the formula to calculate the area from GCS...
    – Pugazh
    Commented Nov 13, 2014 at 6:04
  • I found an answer to your question.
    – xunilk
    Commented Mar 6, 2015 at 8:15
  • I followed steps in Answer above, and I dont have the same result. I have problems with QgsDistanceArea, area is not calculated, or is calculated in degrees. wb=iface.activeLayer() feat = wb.selectedFeatures() area = QgsDistanceArea() #creating object area.setEllipsoid('WGS84') #setting ellipsoid True area.setEllipsoidalMode(True) #setting ellipsoidal mode geom = feat[0].geometry().asPolygon() print "Continental Spain area is %.2f km2" % (area.measurePolygon(geom[0])/1e6) Result is: Continental Spain area is 0.00 km2 I use Qgis Wien 2.8.2
    – Marcel GJS
    Commented Aug 11, 2015 at 7:26

2 Answers 2

4

It's possible with the class QgsDistanceArea of PyQGIS. For testing this, I am going to use the world_borders.shp shapefile and selecting the continental area of Spain (select feature); as in the below image:

enter image description here

After, I ran this code in the Python Console of QGIS:

wb=iface.activeLayer()

feat = wb.selectedFeatures()

geom = feat[0].geometry().asPolygon()

print len(geom[0])  #length of geom (1143 in total)

print geom[0][0:9] #Slice of 9 first QgsPoint (1143 in total)
                   #of polygon geometry (geographic coordinates)

print geom[0][-1] #checking that the last point is the same as
                  #the first one

area = QgsDistanceArea() #creating object

area.setEllipsoid('WGS84')  #setting ellipsoid

area.setEllipsoidalMode(True) #setting ellipsoidal mode

#calculating area
print "Continental Spain area is %.2f km2" % (area.measurePolygon(geom[0])/1e6)

Results were:

Python Console 
Use iface to access QGIS API interface or Type help(iface) for more info
execfile(u'/home/zeito/scriptspyqgis/area_geografic.py'.encode('UTF-8'))
1143  #vertex number for Continental Spain Polygon
[(-0.761016,37.846), (-0.801389,37.8014), (-0.841944,37.7489), (-0.858333,37.7217), 
(-0.86,37.7167), (-0.8575,37.7119), (-0.814167,37.6647), (-0.804722,37.6572), 
(-0.786945,37.6478)]  #Slice of 9 first QgsPoint (1143 in total)
(-0.761016,37.846) #last point
Continental Spain area is 494564.85 km2

Wikipedia reports a total area of 504,782 km2 (including insular territories).

Note editing:

This code takes into account the possibility of Multi parts or Single parts.

area = QgsDistanceArea() #creating object

area.setEllipsoid('WGS84')  #setting ellipsoid

area.setEllipsoidalMode(True) #setting ellipsoidal mode

wb=iface.activeLayer()

feat = wb.selectedFeatures()

multiPart = feat[0].geometry().isMultipart()

if (multiPart is False):
    pol = feat[0].geometry().asPolygon()
    print "area is %.2f km2" % (area.measurePolygon(pol[0])/1e6)
else:
    pol = feat[0].geometry().asMultiPolygon()
    print "There are %d polygons" % len(pol)
    i=1
    for geom in pol:
        print "area%d is %.2f km2" % (i,area.measurePolygon(geom[0])/1e6)
        i +=1

Its implementation in Python Console for Alaska (I used another shapefile) results in:

enter image description here

2
  • any ideea on how to do this in MySQL 5.6 using spatial extensions ? Commented Apr 22, 2015 at 14:00
  • You should call area.computeAreaInit() as mentioned in another answer before calling area.measurePolygon(), otherwise you will get the answer as 0.0 Commented Sep 11, 2015 at 9:19
2

EDITING: CALL QgsDistanceArea().computeAreaInit()


    wb=iface.activeLayer()
    feat = wb.selectedFeatures()
    area = QgsDistanceArea() #creating object
    area.setEllipsoid('WGS84')  #setting ellipsoid
    True
    area.setEllipsoidalMode(True) #setting ellipsoidal mode

    area.computeAreaInit() # CALL THIS BEFORE CALCULATING AREA !!!

    geom = feat[0].geometry().asPolygon()
    print "Continental Spain area is %.2f km2" % (area.measurePolygon(geom[0])/1e6)

Now I have result 494564854638.802 Square Meters (VERY GOOD). Before calculating area, call method of QgsDistanceArea class computeAreaInit().

2
  • 1
    Can you please explain, why you are calling the class computeAreaInit() before calculating the area! This would be interesting for those who want to benefit from your answer.
    – Stefan
    Commented Aug 11, 2015 at 8:59
  • without computeAreaInit, result from Answer for spain feature is 0.00 km2.
    – Marcel GJS
    Commented Aug 12, 2015 at 10:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.