I wrote the following script, which runs fine from the QGIS python console (it loads two shapefiles, intersects them, and then finds the areas of the new regions):

    import qgis.utils
    import qgis.core
    from qgis.analysis import *
    from PyQt4.QtCore import QVariant
    import os
    print(os.getcwd())

    layer1 = iface.addVectorLayer("/Users/casta2k/ShapefilePractice/precincts/USA_precincts.shp", "precincts", "ogr")
    layer2 = iface.addVectorLayer("/Users/casta2k/ShapefilePractice/data_EPSG_4326/Bnd_2015_q1_region.shp", "zipcodes", "ogr")

    overlayAnalyzer = QgsOverlayAnalyzer() 
    overlayAnalyzer.intersection(layer1, layer2, "/Users/casta2k/ShapefilePractice/intersect/intersectOutput.shp")

    layer3 = iface.addVectorLayer("/Users/casta2k/ShapefilePractice/intersect/intersectOutput.shp", "intersect", "ogr")

    layer = iface.activeLayer()
    provider = layer.dataProvider()

    areas = [ feat.geometry().area() 
          for feat in layer.getFeatures() ]

    field = QgsField("area", QVariant.Double)
    provider.addAttributes([field])
    layer.updateFields()

    idx = layer.fieldNameIndex('area')

    for area in areas:
        new_values = {idx : float(area)}
        provider.changeAttributeValues({areas.index(area):new_values})


I'd like to be able to run the code directly from a python script, without having to open the GUI. I've tried searching the internet, but I've found nothing helpful. Do I have to start by running QgsApplication.initQgis()? I thought I had to, but I get a segmentation fault. Also, I tried changing iface.addVectorLayer() to QgsVectorLayer(), which seems to work fine in the Python console, but it seems like the files aren't loaded when I run the python script directly.

As you can see, I'm very lost here (the documentation that is available is very poor), so I would appreciate any help. Thanks!