Note: This anwser is based on gandalf3 correct solution and is related to @skuwid comment in the gandalf3 anwser comments
Here is what I created based on gandalf3 answer to iterate through all selected objects:
import bpy
import math
# Target rotation:
target_rotation = [90, 90, 90]
# Objects to operate on (all selected objects):
selected_objects = bpy.context.selected_objects
# Loop through all selected objects
for ob in selected_objects:
# Deselect all objects first to avoid any overlap
bpy.ops.object.select_all(action='DESELECT')
# Select the current object and make it active
ob.select_set(True)
bpy.context.view_layer.objects.active = ob
# Apply current rotation to zero out before applying new rotation
# Using isolate_users here is optional - in my case I had to use it because children of my objects have been using Multi User Data which dissallowed applying transform
bpy.ops.object.transform_apply(location=False, rotation=True, scale=False, isolate_users = True)
# Get inverse of the rotation we want
t_rot = (math.radians(target_rotation[0] * -1),
math.radians(target_rotation[1] * -1),
math.radians(target_rotation[2] * -1))
# Set rotation mode and apply the inverse rotation in ZYX order
ob.rotation_mode = 'ZYX'
ob.rotation_euler = t_rot # Set rotation to inverse
bpy.ops.object.transform_apply(location=False, rotation=True, scale=False)
# Now apply the actual target rotation
t_rot = (math.radians(target_rotation[0]),
math.radians(target_rotation[1]),
math.radians(target_rotation[2]))
ob.rotation_mode = 'XYZ' # Set rotation mode back to XYZ
ob.rotation_euler = t_rot # Set rotation
# Deselect the current object before moving to the next one
ob.select_set(False)
I'm also not using blender and python very often so beware of any mistakes I may have accidentally put there ;)