I've had the same issue in my work, and ultimately turned to NetCDF utilities to solve it. I know nothing about NetCDF really, just that I can
gdal_translate a raster to that format, and that some of the NetCDF CLI tools can be used to accomplish these "joins".
In my limited testing, the NetCDF approach takes about 40 seconds (most of which is writing the CSV out), while the gdal_translate + sqlite approach takes about 5 and a half minutes.
Here's an example:
import xarray as xr
# convert input tifs to netcdf
subprocess.check_call(['gdal_translate', '-of', 'NetCDF', 'biomass.tif', 'biomass.nc'])
subprocess.check_call(['gdal_translate', '-of', 'NetCDF', 'extent.tif', 'extent.nc'])
# rename bands to match data
subprocess.check_call(['ncrename', '-v', 'Band1,biomass', 'biomass.nc'])
subprocess.check_call(['ncrename', '-v', 'Band1,extent', 'extent.nc'])
# run a netcdf append command to add biomass data to the extent netcdf
subprocess.check_call(['ncks', '-A', 'biomass.nc', 'extent.nc'])
# open extent and then write it directly to CSV
# unfortunately there doesn't seem to be a fast CLI tool for this
ds = xr.open_dataset('extent.nc')
df = ds.to_dataframe().reset_index()
AFAIK there's no way to write one of these "joined" NetCDFs directly to CSV from the command line, but xarray and pandas can read them and write them fine. Hope this helps!