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?

  • 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 Nov 12 '14 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 Nov 12 '14 at 14:45
  • I required the formula to calculate the area from GCS... – Pugazh Nov 13 '14 at 6:04
  • I found an answer to your question. – xunilk Mar 6 '15 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 Aug 11 '15 at 7:26
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

  • any ideea on how to do this in MySQL 5.6 using spatial extensions ? – Stefan Rogin Apr 22 '15 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 – spatialthoughts Sep 11 '15 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().

  • 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 Aug 11 '15 at 8:59
  • without computeAreaInit, result from Answer for spain feature is 0.00 km2. – Marcel GJS Aug 12 '15 at 10:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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