Run my QGIS Python script in headless mode

Question

I have the following Python script that I've been manually pasting into QGIS in order to do some importing/analysis:

import math
import csv
from pathlib import Path
This script ranks properties in Toronto by distance to the nearest rail transit station.
It relies on shapefiles from the following sources:
https://ckan0.cf.opendata.inter.prod-toronto.ca/ne/dataset/ttc-subway-shapefiles
https://ckan0.cf.opendata.inter.prod-toronto.ca/ne/dataset/address-points-municipal-toronto-one-address-repository
class Ranker(object): 
    @staticmethod
    def haversine(lon1, lat1, lon2, lat2):
        """
        Calculate the great circle distance between two points 
        on the earth (specified in decimal degrees)
        """
        # convert decimal degrees to radians 
        lon1, lat1, lon2, lat2 = map(math.radians, [lon1, lat1, lon2, lat2])
        # haversine formula 
        dlon = lon2 - lon1 
        dlat = lat2 - lat1 
        a = math.sin(dlat/2)2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon/2)2
        c = 2 * math.asin(math.sqrt(a)) 
        km = 6367 * c
        return km
@staticmethod
def calculateDistance(featureProperty, layerStationEntrances):
    if featureProperty.geometry() is None:
        print("Null property geometry! Returning -1!")
        return -1
    shortestDistance = 1000
    featuresStationEntrances = layerStationEntrances.getFeatures()
    for featureStationEntrance in featuresStationEntrances:
        if featureStationEntrance.geometry() is None:
            print("Null station entrance geometry! Skipping!")
            continue
        latStationEntrance = featureStationEntrance.geometry().centroid().asPoint().y()
        lonStationEntrance = featureStationEntrance.geometry().centroid().asPoint().x()
        latProperty = featureProperty.geometry().centroid().asPoint().y()
        lonProperty = featureProperty.geometry().centroid().asPoint().x()
        currentDistance = Ranker.haversine(lonStationEntrance, latStationEntrance, lonProperty, latProperty)
        if currentDistance < shortestDistance:
            shortestDistance = currentDistance
    return shortestDistance

@staticmethod
def importAndAnalyze(rowLimit = 3):
    layerProperties = iface.addVectorLayer("/Users/mrmayor/Downloads/municipal-address-points-wgs84-latitude-longitude/ADDRESS_POINT_WGS84.shp", "Toronto_Addresses", "ogr")
    if not layerProperties:
        print("layerProperties failed to load!")
    layerStationEntrances = iface.addVectorLayer("/Users/mrmayor/Downloads/ttc-subway-shapefile-wgs84/TTC_SUBWAY_LINES_WGS84.shp", "Subway_Stations", "ogr")
    if not layerStationEntrances:
        print("layerStationEntrances failed to load!")
    features = layerProperties.getFeatures()
    counter = 0
    featuresSelected = []
    Path('/Users/mrmayor/Desktop/toronto-property-rankings.csv').touch()
    with open('/Users/mrmayor/Desktop/toronto-property-rankings.csv', 'a') as csvFile:
        csvWriter = csv.writer(csvFile, delimiter=',', quotechar='"', quoting=csv.QUOTE_NONNUMERIC)
        csvWriter.writerow(['ID', 'Lat','Lon', 'Distance'])
        for feature in features:
            if counter < rowLimit and feature.geometry().centroid() is not None:
                try:
                    distance = Ranker.calculateDistance(feature, layerStationEntrances)
                    csvWriter.writerow([feature['ID'], feature.geometry().centroid().asPoint().y(),feature.geometry().centroid().asPoint().x(), distance])
                    if counter % 50 == 0:
                        csvFile.flush()
                except ValueError:
                    print("Property could not be converted to point.")
            else:
                break
            counter += 1
        csvFile.flush()


Ranker.importAndAnalyze(10000000)

While this works, in that I get a nice CSV file as output, showing how far each Toronto city property is from its nearest heavy rail transit station, I am a bit annoyed by the need to launch QGIS, create an empty project, paste the Python script into the QGIS console and click the "run" button. Since I'm passing this off to someone else in order to check my work, I don't want a lot of room for ambiguity or human error.

Is there a way -- preferably using a Docker image -- that I can launch this script -- which I have saved in a GitHub repository as rank-and-export.py -- in headless mode from the command line, using easily repeatable steps? If so, how would I go about doing that?