So first, represent your system as the following in your case:
\begin{align}
\frac{\partial q_k}{\partial t} - q_k \nabla \cdot F_k(\boldsymbol{q}) &= 0 & \forall k \in \lbrace 1, 2 \rbrace
\end{align}
where $\boldsymbol{q} = [q_1, q_2]^T = [f, g]^T$, $F_1(\boldsymbol{q}) = q_2 \hat{e}_1$, and $F_2(\boldsymbol{q}) = q_1 \hat{e}_1$. Then you just perform the normal Galerkin for each equation using a set of weight functions for that equation. Let us assume we have some space of functions defined for each function/equation, call them $\mathcal{V}_k$. Then we can proceed with the Galerkin projection like so:
\begin{align}
\int_{\Omega} w_j \left(\frac{\partial q_k}{\partial t} - q_k \nabla \cdot F_k(\boldsymbol{q})\right) d\Omega &= 0 &\forall w_j \in \mathcal{V}_k \\
%
\int_{\Omega} w_j \frac{\partial q_k}{\partial t} d\Omega &= \int_{\Omega} w_j q_k \nabla \cdot F_k(\boldsymbol{q}) d\Omega &\forall w_j \in \mathcal{V}_k \\
%
&= \int_{\Omega} \nabla \cdot \left( w_j q_k F_k(\boldsymbol{q}) \right) - F_k(\boldsymbol{q}) \cdot \nabla \left( w_j q_k \right)d\Omega &\forall w_j \in \mathcal{V}_k \\
%
&= \int_{\Gamma} w_j q_k F_k(\boldsymbol{q}) \cdot \hat{n}\; d\Gamma - \int_{\Omega} F_k(\boldsymbol{q}) \cdot \nabla \left( w_j q_k \right)d\Omega &\forall w_j \in \mathcal{V}_k \\
&= G_k(\boldsymbol{q}, w_j) &\forall w_j \in \mathcal{V}_k
\end{align}
From here you can discretize things further with your approximate representation for the functions $q_k(x)$ using your space of functions $\mathcal{V}_k$ and come up with the equations you need in terms of any coefficients you are solving for.