How to get local rotation from Matrix but taking parent transforms into account
So quaternions worked! Thanks BatFinger!

import bpy
from mathutils import Vector

Original_rot = bpy.data.objects['Armature'].pose.bones['shoulder'].matrix_basis.copy()

bpy.data.objects['Armature'].pose.bones['shoulder'].rotation_euler = Vector((0,0,0))
bpy.context.evaluated_depsgraph_get()


a = bpy.data.objects['Armature'].pose.bones['shoulder'].matrix.copy()
b = bpy.data.objects['Armature'].pose.bones['shoulder_extra_snap'].matrix.copy()

bpy.data.objects['Armature'].pose.bones['shoulder'].matrix = b

q1 = a.to_quaternion()

q2 = b.to_quaternion()

difQuat = q1.rotation_difference(q2)

difEuler = difQuat.to_euler()

bpy.data.objects['Armature'].pose.bones['shoulder'].rotation_euler = Vector(difEuler) + Vector(Original_rot.to_euler())


Might not be the best code but it works!