The reason you don't get conservation is that you've **used the product rule before discretizing**, so conservation would require an exact cancellation of truncation errors in the different product terms (which generally won't happen).  

Instead, you should **directly discretize the conservative form of the equation**.  Using the usual notation for finite differences and setting $(\gamma^x,\gamma^y)=(u,v)$, this would be something like:

$$\theta^{n+1}_{i,j} = \theta^{n-1}_{i,j} - \frac{\Delta t}{\Delta x}(\theta^n_{i+1,j}u^n_{i+1,j} - \theta^n_{i-1,j}u^n_{i-1,j}) - \frac{\Delta t}{\Delta y}(\theta^n_{i,j+1}v^n_{i-1,j+1} - \theta^n_{i,j}u^n_{i-1,j}).$$

It's easy to check that if you sum $\theta^{n+1}$ over the whole grid, all of the fluxes cancel out except for those at the boundaries.  Here I've used centered differences in time and space.  This will be fine if your initial data is smooth and well-resolved, and if second-order accuracy (in $\Delta t$ and $\Delta x$) is sufficient.  For more complicated or demanding situations, you might refer to (for instance) [Randall LeVeque's *Finite Volume Methods for Hyperbolic Problems*](http://depts.washington.edu/clawpack/book.html) and the [Clawpack software](http://www.clawpack.org).

Note: In the literature I'm familiar with, the conservation law you've written **is** referred to as the 2D advection equation.  The non-conservative version is sometimes given the same name, or may be referred to as the **color equation**.