Is there a python option to "join attributes by location"?

You may want to take a look at Shapely and Fiona. Fiona is a wrapper for gdal to make spatial file import and export easy. Shapely provides geometry functionality. Here is a very simple example to give you the idea. It joins polygon attributes to all points within that polygon.

The example data I have used are these polygons and these points.

import fiona
from shapely.geometry import shape
from copy import deepcopy

with"planning_neighborhoods.shp", "r") as n: 

    with"Schools_Private_Pt.shp", "r") as s:

        # create a schema for the attributes
        outSchema =  deepcopy(s.schema)

        with ("Schools_withNbhd.shp", "w", s.driver, outSchema, as output:

            for school in s: 
                for neighborhood in n:
                    # check if point is in polygon and set attribute
                    if shape(school['geometry']).within(shape(neighborhood['geometry'])):  
                        school['properties']['neighborho'] = neighborhood['properties']['neighborho'] 
                    # write out
                            'properties': school['properties'], 
                            'geometry': school['geometry']

Although still a bit rough around the edges, especially when it comes to documentation and examples, but geopandas' future looks bright. It basically combines the power of pandas dataframes with geospatial capabilities of shapely.

the function you look for is called sjoin

Make sure your machine/instance has enough memory to perform the operation

import geopandas as gpd
import pandas as pd
import os

gdfLeft = gpd.read_file(os.path.join(PATH,INPUT_FILE_NAME_1))
gdfRight = gpd.read_file(os.path.join(PATH,INPUT_FILE_NAME_2))

gdfJoined = gpd.sjoin(gdfLeft, gdfRight, how="left", op='intersects')