I do not know if this will be helpful at all but this is how I am doing my masks and population:
population_image_slice_arr, new_image_transform = rasterio.mask.mask(world_pop_image, [polygon_for_some_country], crop=True, nodata=0)
population_image_slice_arr[population_image_slice_arr < 0] = 0
total_population = population_image_slice_arr.sum()
I know this is incorrect because when i did manual scaling using the bounds of the image and shapely.ops.transform the population for say India would be 800 million when it should be 1.3 billion but then when I changed it to use the affine matrix the math came out to 0.0 population.