You wouldn't need to convert all the shapefiles to raster at once, so volume consideration is a non-issue.
Using a raster masking technique is going to result in imprecision at the pixel boundary -- Does a file which include the center of the pixel count as "covered" if not 100% of the pixel is covered? Using a vector technique could attain higher precision, at the cost of a lot of processing.
In raster space, you just need to mask each vector at a time (with consideration as to whether you should take a one pixel (or
sqrt(2)*pixeldim) negative buffer before masking), sum the presence/absence flag with a cumulator, and then divide by the number of shapefiles (encoding this as byte (either 0-100 or 0-250) would speed future processing, at the cost of precision)
In vector space, you'd need to:
- Determine the resolution of your percentage raster (pixel size, type and depth)
- Create a partition vector to form the basis for subsetting the vector space which is in alignment with the raster (large enough to reduce the iterations, but small enough to save on vector overlay cost)
- Iterate over the tiles with a clip and union of each shapefile, then compute percentage and convert that to raster, encode the desired pixel value, and mask this into a final raster.
For most approaches to either technique, you're going to need a full GIS, not just GDAL and Python.