Calculate area in square meters from degrees?

Question

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?

Answer 1

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:

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:

Answer 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().