Lat / Lon are in degrees
    Depth in km
    quantity_ID needs to match those in the litho1 model

    Points that are not found are given the out-of-range value of -99999

def property_at_lat_lon_depth_points(lat, lon, depth, quantity_ID="DENSITY"):
    """
    Lat / Lon are in degrees
    Depth in km
    quantity_ID needs to match those in the litho1 model

    Points that are not found are given the out-of-range value of -99999
    """

    nlayers = len(l1_layer_decode)
    shape = np.array(lon).shape

    lon1   = np.array(lon).reshape(-1)
    lat1   = np.array(lat).reshape(-1)
    depth1 = np.array(depth).reshape(-1)
    point_properties = np.empty_like(depth1)

    layer_depths     = np.empty((nlayers, lat1.shape[0]))
    layer_properties = np.ones((nlayers+1, lat1.shape[0])) * -99999.0  # if the point is not found it will end up in the overshoot !

    # should assert here that the three arrays are equal size

    for i in range(0, nlayers, 1 ):
        layer_depths[i], err = _interpolator.interpolate( lon1 * np.pi / 180.0, lat1 * np.pi / 180.0,
                                      _litho_data[i,l1_data_decode["DEPTH"]], order=1)
        layer_properties[i], err = _interpolator.interpolate( lon1 * np.pi / 180.0, lat1 * np.pi / 180.0,
                                      _litho_data[i,l1_data_decode[quantity_ID]], order=1)


    A = -layer_depths
    B = -depth1 * 1000.0
    C = divmod(np.searchsorted(A.ravel(), B), A.shape[1])[0] # YEP - this seems to be the best way !!

    # point_properties = np.diag(layer_properties[C[:],:])

    point_properties = np.empty_like(depth1)
    for i,layer in enumerate(C):

        point_properties[i] = layer_properties[layer,i]

    return C, point_properties.reshape(shape)