How to find out whether a point coordinate exists on a path between other two point coordinates?

Question

I've created an application that scrapes data on house listings and calculates the distance to the nearest shop and train station. I've got the lat,lon coordinates for all the shops, listings and train stations and I currently use geopy.distance

import geopy.distance
distanceToShop = geopy.distance.distance(shopCoordinate, listingCoordinate).km
distanceToStation = geopy.distance.distance(stationCoordinate, listingCoordinate).km

These yield a value in kilometers between the two points. Then I filter out locations that are not close enough to a shop and trainstation. For train stations I tolerate a longer distance than for shops. Sometimes the shop is too far to be filtered in, but on the way to the station. In situations like this I would tolerate a longer distance to the shop, since it's on the way to the station.

In this case the shop is between the house and the stations. Station is at coordinate (40.651100, -73.949626). House is at coordinate (40.651401, -73.937987). Distance between house and station is 0.98 kilometers.

In the below scenario shop is at coordinate (40.650913, -73.948301) and distance from house to shop is 870 meters.

Acceptable distance from house to shop is 650 meters. However, if shop is on the way from the station a longer distance is acceptable. If I naively drop out houses because the closest shop is too far (870m > 650m) then this house would not qualify, but in this case the shop is on the way from the station (on the road home), so it should qualify.

Representation of problem on a map

A more robust solution would calculate an area between the two points (e.g. a square), and find out whether the shop is located within the square. This would allow for a more realistic result.

Representation of square drawn on map

How can I find out programatically if this is the case? Are there libraries that I could study?

EDIT: SOLVED - found out a way to calculate the bearing (degree angle between two coordinates) which can be manipulated and used to calculate other coordinates. By calculating the bearing I am able to offset it by +/- 90 degrees to create corners for both the target and destination coordinates. Then using another import (shapely) I drew a polygon from the coordinates and checked whether the shop existed in that matrix.

EDIT 2: Neater solution using geographiclib.geodesic with fewer imports:

import geopy.distance
from geographiclib.geodesic import Geodesic
from shapely.geometry import Point
from shapely.geometry.polygon import Polygon

def getAreaCorners(origin, target, distanceInMeters, angle):
    geod = Geodesic.WGS84
    bearing = geod.Inverse(origin[0], origin[1], target[0], target[1])['azi1']
    rightBearing = bearing + angle
    if rightBearing < 0: rightBearing += 360
    leftBearing = bearing - angle
    if leftBearing < 0: leftBearing += 360
        
    og1 = geod.Direct(origin[0], origin[1], rightBearing, distanceInMeters)
    og2 = geod.Direct(origin[0], origin[1], leftBearing, distanceInMeters)
    og1 = (og1['lat2'], og1['lon2'])
    og2 = (og2['lat2'], og2['lon2'])
    
    tg1 = geod.Direct(target[0], target[1], rightBearing, distanceInMeters)
    tg2 = geod.Direct(target[0], target[1], leftBearing, distanceInMeters)
    tg1 = (tg1['lat2'], tg1['lon2'])
    tg2 = (tg2['lat2'], tg2['lon2'])
    return [og1, og2, tg2, tg1]

def isShopOnPath(areaCoordinateList, shopCoordinate):
    polygon = Polygon(areaCoordinateList)
    shopCoordinatePoint = Point(shopCoordinate)
    return polygon.contains(shopCoordinatePoint)

origin = (40.650809, -73.949583)
target = (40.644661, -73.943108)
shopCoordinate = (40.650467, -73.948553)
distanceInMeters = 100
angle = 90

distanceToShop = round(geopy.distance.distance(shopCoordinate, target).km,2)
distanceToStation = round(geopy.distance.distance(origin, target).km,2)

print('House is {} km from station. Shop is {} km from house.'.format(distanceToStation, distanceToShop))
if distanceToShop >= 0.65:
    if isShopOnPath(getAreaCorners(origin, target, distanceInMeters, angle), shopCoordinate):
        print('Shop is too far, but on the way from the station.')
    else:
        print('Shop is too far and not on the way from the station')
else:
    print('Shop is less than 650m from the house.')