How to do quantum mechanics if the Schrödinger equation was a Sturm-Liouville equation?
Nowhere in the actual formalism of QM is it necessary that the inner product be the familiar unweighted one given by
$$\langle f,g\rangle = \int_a^b \overline{f(x)}g(x) dx$$
In particular, the definitions of Hermiticity and self-adjointness remain the same. If the system is in some state $\psi$, the probability of measuring some observable $A$ to be in some (Borel) set $E\subseteq \mathbb R$ is given by
$$\operatorname{Prob}_A(E,\psi) = \frac{\Vert \Pi_A(E) \psi\Vert^2}{\Vert \psi\Vert^2}$$
where $\Pi_A(E)$ is the projection onto the relevant spectral subspace and $\Vert \psi\Vert^2 = \langle \psi,\psi\rangle$.
For example, if we're thinking about the position operator $X$, then the projector $\Pi_X(E)$ eats a wavefunction $\psi$ and spits out $\chi_E\cdot \psi$, where $\chi_E(x) = \cases{ 1 & $x\in E$\\0 & else}$ is the indicator function on $E$. As a result, the probability that the particle is in e.g. the interval $[0,1]$ is
$$\operatorname{Prob}_X([0,1],\psi) = \frac{\Vert \chi_{[0,1]}\psi\Vert^2}{\Vert \psi\Vert^2}= \int_0^1 w(x) \overline{\psi(x)}\psi(x) dx$$ assuming that $\psi$ is normalized.
If instead we are interested in an operator $L$ with a discrete spectrum, then it will have an orthonormal eigenbasis $\{\phi_n\}$ where orthonormality is once again defined via the weighted inner product, and so everything proceeds as usual from there.
To directly address your questions,
Can one write the equation in operator form? $$\hat L |\psi_n\rangle = E_n|\psi_n\rangle$$
Yes, this has nothing to do with the inner product.
Is $\mathbb I = \sum_n |\psi_n\rangle\langle \psi_n|$ still true?
Since the eigenstates of a self-adjoint operator $|\psi_n\rangle$ can be orthonormalized (with respect to the weighted inner product), then this remains true.
Related to the above what would $\langle r'|\hat L|r\rangle$ and $\langle r|\psi\rangle$ be? Is $\langle r|r'\rangle=\delta(r-r')$ still true?
This ultimately comes down to how we'd like to scale $|x\rangle$. Since these "states" cannot actually be normalized, we have to choose the form we'd like the so-called resolution of the identity to take.
Convention 1: $\mathbb I = \int |x\rangle\langle x| dx$
Since $\mathbb I^2 = \mathbb I$, this convention would dictate that $\langle x|x'\rangle = \delta(x-x')$, since
$$\mathbb I^2 = \int dx \int dx' |x\rangle\langle x|x'\rangle\langle x'|$$ From there, we would have $$|\psi\rangle = \int |x\rangle \langle x|\psi\rangle dx$$ $$\implies \Vert \psi\Vert^2 = \langle \psi|\psi\rangle = \int dx\int dx' \ \langle \psi|x'\rangle\langle x|\psi\rangle\langle x'|x\rangle$$ $$ = \int dx |\langle x|\psi\rangle |^2 = \int w(x)\overline{\psi(x)}\psi(x) \ dx$$
$$\implies \langle x|\psi\rangle = \sqrt{w(x)} \psi(x)$$
Furthermore, the action of $L$ in the position basis would be
$$\int dx\int dx' |x\rangle\langle x|L|x'\rangle\langle x'|\psi\rangle $$ $$\implies \langle x'|L|x\rangle = -\frac{1}{w(x)}\left( \frac{d^2}{dx^2} - V(x)\right)\frac{1}{\sqrt{w(x)}}\delta(x-x')$$
which is chosen to give the result $$ [\ldots] = \int dx | x\rangle \frac{-1}{w(x)}\left(\frac{d^2}{dx^2} - V(x)\right) \psi(x)$$
Convention 2: $\mathbb I = \int dx\ w(x) |x\rangle\langle x|$
Enforcing $\mathbb I^2 = \mathbb I$ yields the result that $\langle x'|x\rangle = \frac{1}{w(x)}\delta(x-x')$. Inserting the identity operator into $\langle \psi|\psi\rangle$ then gives that $\langle x|\psi\rangle = \psi(x)$, as usual. Lastly, requiring that
$$\int dx\int dx' \ w(x)|x\rangle\langle x|L|x'\rangle w(x')\langle x'|\psi(x)\rangle = \int dx |x\rangle \frac{-1}{w(x)}\left(\frac{d^2}{dx^2} - V(x)\right) \psi(x)$$
yields the result
$$\langle x|L|x'\rangle = -\frac{1}{w(x)^2} \left(\frac{d^2}{dx^2}-V(x)\right) \frac{1}{w(x)}\delta(x-x')$$
You could imagine different conventions, of course. I'd also like to emphasize again that none of this is physically meaningful, in the sense that it all comes down to how you'd like to normalize $|x\rangle$ to make your life most convenient.