Faster way to retrieve subsets of data from large VRT dataset by location/bounding box

Question

What I need:

I need to speed up queries as much as possible to a series of large VRT datasets (total is around 1TB of data, although each VRT is much smaller than that).

I am open to ANY suggestions, including other programming languages, and/or other libraries, or else just improvements to my code below. Possibly what I need is to do some transformations on the VRTs themselves?

Context:

Previously posted this on GIS stack exchange but for some reason they didn't think it fit there, even though it is GIS specific.

Currently I'm using Python with gdal to query these VRTs, here is the key piece of code in question:

# Retrieve a sub-image as a matrix of numbers
# 'boundingBox' is a simple list, coordinates of the corners of the sub-image I would like

def queryVRT(boundingBox, pathToVRT):
        res = 250
        igh = "+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs"
        kwargs = {'format': 'MEM', 'projWin': boundingBox,
                  'projWinSRS': igh, 'xRes': res, 'yRes': res}
        memory_ds = gdal.Translate('', pathToVRT, **kwargs)
        my_array = memory_ds.ReadAsArray()
        del memory_ds
        return my_array

The VRTs in question are actually publicly available here: https://www.isric.org/explore/soilgrids

But here is a small sample of one of the .vrt files, incase it helps (I'm not great at reading them):

<VRTDataset rasterXSize="159243" rasterYSize="58034">
  <SRS>PROJCS["Interrupted_Goode_Homolosine",GEOGCS["GCS_WGS_1984 ellipse",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Interrupted_Goode_Homolosine"],UNIT["Meter",1]]</SRS>
  <GeoTransform> -1.9949000000000000e+07,  2.5000000000000000e+02,  0.0000000000000000e+00,  8.3610000000000000e+06,  0.0000000000000000e+00, -2.5000000000000000e+02</GeoTransform>
  <Metadata>
    <MDI key="Code_version">v2.0.0</MDI>
    <MDI key="Covariates">depth,clm_wcl_s07rad,clm_wcl_srdyrsum,clm_wcl_s08rad,luc_lds_mng,mor_mrg_vdp,clm_wcl_s04rad,clm_wcl_bio16,luc_gfc_trely10,clm_wcl_bio07,mor_env_demm,mor_mrg_twi,clm_esa_s09avg,clm_wcl_srdyrstd,luc_esa_l170,clm_mod_lstd01std,clm_esa_s05avg,clm_wcl_p05tot,mor_mrg_vbf,clm_wcl_p04tot,clm_wcl_bio18,clm_wcl_p10tot,clm_wcl_p09tot,clm_mod_lstd03std,veg_mod_evirng,veg_mod_eviyravg,clm_wcl_bio19,clm_mod_lstd09std,clm_wcl_p12tot,clm_wcl_bio08,clm_wcl_bio17,clm_mod_lstd11std,clm_mod_lstd05std,clm_mod_lstd02std,clm_esa_s08avg,clm_wcl_bio15,clm_mod_lstd10std,veg_mod_nppy15,clm_mod_lstd04std,clm_mod_lstd12std</MDI>
    <MDI key="Depth_interval"></MDI>
    <MDI key="Litter_layers">FALSE</MDI>
    <MDI key="Mapped_units"></MDI>
    <MDI key="Model">Quantile Regression Forests</MDI>
    <MDI key="Model_type">ranger</MDI>
    <MDI key="Mtry">12</MDI>
    <MDI key="Number_trees">200</MDI>
    <MDI key="Other_synth_profiles"></MDI>
    <MDI key="Outputs_abstract">no WISE, March 2020 snapshot, NZ and Siberia gaps</MDI>
    <MDI key="Outputs_version">RUN03</MDI>
    <MDI key="Quantile"></MDI>
    <MDI key="Transformation"></MDI>
    <MDI key="WoSIS_version">Data stream 7</MDI>
  </Metadata>
  <VRTRasterBand dataType="Int16" band="1">
    <NoDataValue>-32768</NoDataValue>
    <ColorInterp>Gray</ColorInterp>
    <ComplexSource>
      <SourceFilename relativeToVRT="1">./bdod_0-5cm_mean/tileSG-018-027/tileSG-018-027_1-1.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="600" RasterYSize="600" DataType="Int16" BlockXSize="600" BlockYSize="6" />
      <SrcRect xOff="0" yOff="0" xSize="600" ySize="600" />
      <DstRect xOff="48246" yOff="31441" xSize="600" ySize="600" />
      <NODATA>-32768</NODATA>
    </ComplexSource>
    <ComplexSource>
      <SourceFilename relativeToVRT="1">./bdod_0-5cm_mean/tileSG-018-027/tileSG-018-027_3-3.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="600" RasterYSize="600" DataType="Int16" BlockXSize="600" BlockYSize="6" />
      <SrcRect xOff="0" yOff="0" xSize="600" ySize="600" />
      <DstRect xOff="49446" yOff="32641" xSize="600" ySize="600" />
      <NODATA>-32768</NODATA>
    </ComplexSource>
    <ComplexSource>
      <SourceFilename relativeToVRT="1">./bdod_0-5cm_mean/tileSG-018-027/tileSG-018-027_3-1.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="600" RasterYSize="600" DataType="Int16" BlockXSize="600" BlockYSize="6" />
      <SrcRect xOff="0" yOff="0" xSize="600" ySize="600" />
      <DstRect xOff="48246" yOff="32641" xSize="600" ySize="600" />
      <NODATA>-32768</NODATA>
    </ComplexSource>
    <ComplexSource>
      <SourceFilename relativeToVRT="1">./bdod_0-5cm_mean/tileSG-018-027/tileSG-018-027_3-2.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="600" RasterYSize="600" DataType="Int16" BlockXSize="600" BlockYSize="6" />
      <SrcRect xOff="0" yOff="0" xSize="600" ySize="600" />
      <DstRect xOff="48846" yOff="32641" xSize="600" ySize="600" />
      <NODATA>-32768</NODATA>
    </ComplexSource>
    <ComplexSource>
      <SourceFilename relativeToVRT="1">./bdod_0-5cm_mean/tileSG-018-027/tileSG-018-027_2-3.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="600" RasterYSize="600" DataType="Int16" BlockXSize="600" BlockYSize="6" />
      <SrcRect xOff="0" yOff="0" xSize="600" ySize="600" />
      <DstRect xOff="49446" yOff="32041" xSize="600" ySize="600" />
      <NODATA>-32768</NODATA>
    </ComplexSource>

.... etc