dem = rasterio.open('merged.vrt')
g = gpd.read_file('domains_tile_shapefile/hou_grid.geojson').to_crs(dem.crs)

def zonal_average(geo):
    try:
        arr, aff = mask(dem, [geo.geometry], crop=True, all_touched=True, nodata=dem.nodata)
        arr = arr[arr>-9000]
        return np.mean(arr)
    except:
        return None

g['dem_avg'] = g.apply(zonal_average, axis=1)