Thanks Alexey. You are the best! I definitely would not have been able to work that out myself.
For other readers of the forum: here's my final working code:
# Define: Map region
#MapRegionVec = [CenterX, CenterY, SizeX0, SizeY0, ZRotationDeg]
MapRegionVec = [304331.4, 5740607.3, 173.8, 143.6, -25.0]
# ALIGN PHOTOS
chunk.alignPhotos()
# SET MAPPING REGION (BOUNDING BOX)
# create new region object
newregion = PhotoScan.Region()
# Set region center:
centerUTM = PhotoScan.Vector([MapRegionVec[0], MapRegionVec[1], 0])
centerGEO = chunk.crs.unproject(centerUTM)
centerGEO.size = 4
centerGEO.w = 1
centerLocal = chunk.transform.inv() * centerGEO
newregion.center = PhotoScan.Vector([centerLocal[0], centerLocal[1], chunk.region.center[2]])
# Set region size
newregion.size = PhotoScan.Vector([MapRegionVec[2], MapRegionVec[3], chunk.region.size[2]])
# Set region rotation
# build rotation matrix in our coordinate system
RotZDeg = -MapRegionVec[4]
SinRotZ= math.sin(math.radians(RotZDeg))
CosRotZ= math.cos(math.radians(RotZDeg))
RotMat = PhotoScan.Matrix([[CosRotZ,-SinRotZ,0,0],[SinRotZ,CosRotZ,0,0],[0,0,1,0],[0,0,0,1]])
# rotate region bounding box
T = chunk.transform
v = PhotoScan.Vector( [0,0,0,1] )
v_t = T * v
v_t.size = 3
m = chunk.crs.localframe(v_t)
m = m * T
m = RotMat*m
s = math.sqrt(m[0,0]**2 + m[0,1]**2 + m[0,2]**2) #scale factor
R = PhotoScan.Matrix( [[m[0,0],m[0,1],m[0,2]], [m[1,0],m[1,1],m[1,2]], [m[2,0],m[2,1],m[2,2]]])
R = R * (1. / s)
newregion.rot = R.t()
# put newregion to chunk
chunk.region = newregion