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