This is the expected value of a function of the Normal random variable $X \sim {\rm N}(0, \sigma)$, which more over is truncated from below. By the definition then,
\begin{align}
I_1= E\left[\Phi\left(\frac{aX + b}{c}\right) | X > k\right] &= \frac{\int_{k}^{\infty}\frac 1 {\sigma}\phi\left(x/\sigma\right)\Phi\left(\frac{ax + b}{c}\right)dx}{1-\Phi(k/\sigma)}\\
&\\
&\implies I_1 \equiv \frac{I_2}{1-\Phi(k/\sigma)}. \tag{1}
\end{align}
2. We consume the scale factor $\sigma$, by setting $z=x/\sigma$, to get
$$I_2 = \int_{k/\sigma}^{\infty}\phi\left(z\right)\Phi\left(\frac{a\sigma}{c}z + \frac bc\right)dz. \tag{2}$$
3. We express the definite integral $I_2$ as a difference,
\begin{align}
I_2 \equiv I_3-I_4 &\implies \int_{k/\sigma}^{\infty}\phi\left(z\right)\Phi\left(\frac{a\sigma}{c}z + \frac bc\right)dz \\
&= \int_{-\infty}^{\infty}\phi\left(z\right)\Phi\left(\frac{a\sigma}{c}z + \frac bc\right)dz - \int_{-\infty}^{k/\sigma}\phi\left(z\right)\Phi\left(\frac{a\sigma}{c}z + \frac bc\right)dz \tag{3}
\end{align}
4. For integral $I_3$ we have from Owen, D. B. (1980). A table of normal integrals: Communications in Statistics-Simulation and Computation, 9(4), 389-419., p. 403, eq. 10,010.8
$$I_3=\int_{-\infty}^{\infty}\phi\left(z\right)\Phi\left(\frac{a\sigma}{c}z + \frac bc\right)dz = \Phi\left(\frac{b}{\sqrt{c^2+a^2\sigma^2}}\right). \tag{4}$$
5. We turn to $I_4$. From Azzalini, A. and Ant. Capitanio (2014). The skew-normal and related families Cambridge University Press. we have the definition of the "Extended Skew Normal" distribution (p. 36 eq. 2.39) that has density
$$f(z) = \phi(z)\frac{\Phi\left(\gamma z + \tau\sqrt{1+\gamma^2}\right)}{\Phi(\tau)} \tag{5}$$
By definition its CDF, say for quantile $k/\sigma$ is
$$F_Z(k/\sigma) \equiv \Pr(Z \leq k/\sigma) = \int_{-\infty}^{k/\sigma} \phi(z)\frac{\Phi\left(\gamma z + \tau\sqrt{1+\gamma^2}\right)}{\Phi(\tau)}dz.\tag{6}$$
Azzalini & Capitanio (p. 40 eq. 2.48) give this CDF in terms of the Bivariate standard Normal integral $\Phi_2$ as
$$\int_{-\infty}^{k/\sigma} \phi(z)\frac{\Phi\left(\gamma z + \tau\sqrt{1+\gamma^2}\right)}{\Phi(\tau)}dz = \frac{\Phi_2\left(k/\sigma,\,\tau\,;\, -\delta\right)}{\Phi(\tau)},\quad \delta = \frac{\gamma}{\sqrt{1+\gamma^2}}. \tag{7}$$
Here $-\delta$ is the correlation coefficient of the bivariate Normal vector. From this we get
$$\int_{-\infty}^{k/\sigma} \phi(z)\Phi\left(\gamma z + \tau\sqrt{1+\gamma^2}\right)dz = \Phi_2\left(k/\sigma,\,\tau\,;\, -\delta\right). \tag{8}$$
Combining all results we arrive at, for $X \sim {\rm N}(0,\sigma)$,
$$E\left[\Phi\left(\frac{aX + b}{c}\right) | X > k\right] = \frac{\Phi\left(\frac{b}{\sqrt{c^2+a^2\sigma^2}}\right) - \Phi_2\left(k/\sigma,\,\tau\,;\, -\delta\right)}{1-\Phi(k/\sigma)}$$
with
$$\gamma = \frac{a\sigma}{c},\quad \tau\sqrt{1+\gamma^2} = \frac bc ,\quad \delta = \frac{\gamma}{\sqrt{1+\gamma^2}}.$$
6. Mapping the auxiliary coefficients to the original vector $(a,b,c,k)$ of parameters, we finally get
$$E\left[\Phi\left(\frac{aX + b}{c}\right) | X > k\right] = \frac{\Phi\left(\frac{b}{\sqrt{c^2+a^2\sigma^2}}\right)-\Phi_2\left(k/\sigma,\,\frac{b}{\sqrt{c^2+a^2\sigma^2}}\,;\, -\delta\right)}{\Phi(-k/\sigma)},$$
$$\delta = \frac{a\sigma}{\sqrt{c^2+a^2\sigma^2}}. \tag{9}$$
Indicative simulations verify the above result.
Other integrals of similar form can be found in https://en.wikipedia.org/wiki/List_of_integrals_of_Gaussian_functions