Glad you are making progress.
You should be able to use the GUI to create “mask” bands from the shapefile that can later be used with band-maths to set pixels outside the region to a “no data” value. I should also mention that if you end goal is areal statistics for the regions, then you might be to use StatisticsOP.
I wonder if some of what you are doing is more in the domain of a GIS package such as QGIS or ArcGIS (both support Python scripts).