Inconsistently getting "Input shapes do not overlap raster", CRS appear to match

Masking topographic raster data by a shapefile (bedrock geology polygon).

Following this guide I do this:

with fiona.open("geology.shp", "r") as shapefile:
     geology_mask = [feature["geometry"] for feature in shapefile]
out_image, out_transform = rasterio.mask.mask(raster, geology_mask, crop=True)

And it had been working so well on the rasters I plucked from my state DEMs. Now I’ve grabbed many tiles of DEMs and want to iterate through each one to collect data masked by the shapefile. I’m finding the code works (perfectly well!) for a single one of those DEMs; the other 35 give me the error:

     ValueError: Input shapes do not overlap raster.

I have read many answers in which the coordinate reference systems might be off which could explain why they’d overlap in Arc but not for the script. I checked the metadata. Fiona says this about my shapefile:

  "crs": "EPSG:26918",
  "crs_wkt": "PROJCS["NAD83 / UTM zone 18N",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-75],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","26918"]]",

And my TIFFs have the same exact CRS and projections (I extracted them from the same big DEM!) – which is why it’s confusing for me that I might have some CRS issue if everything is the same!:

PROJCRS["NAD83 / UTM zone 18N",
BASEGEOGCRS["NAD83",
    DATUM["North American Datum 1983",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4269]],
CONVERSION["UTM zone 18N",
    METHOD["Transverse Mercator",
        ID["EPSG",9807]],
    PARAMETER["Latitude of natural origin",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8801]],
    PARAMETER["Longitude of natural origin",-75,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8802]],
    PARAMETER["Scale factor at natural origin",0.9996,
        SCALEUNIT["unity",1],
        ID["EPSG",8805]],
    PARAMETER["False easting",500000,
        LENGTHUNIT["metre",1],
        ID["EPSG",8806]],
    PARAMETER["False northing",0,
        LENGTHUNIT["metre",1],
        ID["EPSG",8807]]],
CS[Cartesian,2],
    AXIS["(E)",east,
        ORDER[1],
        LENGTHUNIT["metre",1]],
    AXIS["(N)",north,
        ORDER[2],
        LENGTHUNIT["metre",1]],
USAGE[
    SCOPE["unknown"],
    AREA["North America - 78°W to 72°W and NAD83 by country"],
    BBOX[28.28,-78,84,-72]],
ID["EPSG",26918]]

I have tried to tests, both to no avail:

1. I tried to transform one tile’s srs to match the "working" tile’s srs

   gdalsrsinfo -o wkt goodtile.tif > target.wkt

followed by a gdalwarp

2. I saw a similar question where the CRS were different between the raster and the shapefile and while I don’t think that was the case, I tried to change the CRS for my raster to match my shapefile (like I said I think everyone is in EPSG:26918 already?)

   gdalwarp -t_srs EPSG:26918 badtile.tif badtile_proj.tif

I’m really bothered by the fact that it works for ONE DEM (and others I’d been working with) but not others that are otherwise exactly the same. Any insight?