How to interpolate values in non-rectangular coordinates using Python?

I need to make compensation of values in my optical system using Python. I've measured the dependency of my compensation params in the corners of my table and I want to interpolate such value there linearly, but a map is not a rectangle.
Example:

# Corners coordinates:
a_real = (45, 45) 
a_coeff = (333, 223)

b_real = (-45, -45)
b_coeff = (325, 243)

c_real = (-45, 45)
c_coeff = (339, 244)

d_real = (45, -45)
d_coeff = (319, 228)

Let's say, I want to know compensation coefficients in points (40, 40), or (0, 0).

1. How this can be done? I'm looking at scipy.interpolate.interp2d but I'm not sure that it is my case
2. What if I want to add more points, defining my grid?

You are dealing with 2D unstructured points, so you could use, for instance, the SciPy's interpolate.griddata function.

I would suggest the following solution. I have rearranged your data in a more convenient way. You can add all the points you want to the points NumPy array.

import numpy as np
from scipy.interpolate import griddata

points = np.array([[333, 223], [325, 243], [339, 244], [319, 228]])

values_a = np.array([45, -45, -45, 45])
values_b = np.array([45, -45, 45, -45])
new_value_a = griddata(points, values_a, (325, 232), method="nearest")
new_value_b = griddata(points, values_b, (325, 232), method="nearest")

The code allows you to compute the two values for a new point having coordinates (325, 232). If the new point lies outside the convex hull defined by your points, then you need to set the value of the fill_value parameter (unless you are using the nearest method as explained in the documentation).