What I have a hard time understanding is the connection between line integrals of vector fields and Greens theorem.
It was explained that taking lines integrals of parametrized curves is to be interpreted as the work done over the curve.
However what Greens theorem provides is that we could write this in-terms of a double integral of the area. Bounded by a closed curve.
What are we then calculating? Is it still the work over a closed curve?
As Thom says in the comments, one of the better interpretations is in terms of electrical potential. Suppose that ${\bf E}$ is the electric field due to a charge in a region $R$. One might try to measure the flux of this field through the boundary, denoted $\partial R$ of the region $R$. What this means is that for each point in the boundary, one computes the amount of ${\bf E}$ which flows through this curve and then sums along the curve. If ${\bf n}$ is a unit outward normal vector field along $\partial R$, then this amount is ${\bf E} \cdot {\bf n}$, a scalar quantity that can be summed along the curve, i.e. the line integral $\oint_{\partial R} {\bf E} \cdot {\bf n} \, ds,$ where $ds$ is the element of arc length.
On the other hand, the divergence of ${\bf E}$ at a point can be interpreted as the amount of outward flow through a small circle around this point. If the divergence is negative, then the flow is inward. (This is a little out of order logically, since Green's theorem typically is used by mathematicians to justify this interpretation of the divergence. But my impression is that physicists take this interpretation as almost an intuitive definition of the divergence.) If we sum the divergence over all the points of the region $R$, then we obtain the double integral $\iint_R {\rm div} \, {\bf E} \, dA$, where $dA$ is the element of area.
That these two integrals are equal is the content of Green's theorem. In calculus books, the equation in Green's theorem is often expressed as follows:
$$\oint_{C} {\bf F} \cdot d{\bf r} = \iint_{R} \left( \frac{\partial N}{\partial x} - \frac{\partial M}{\partial y} \right) \, dA,$$
where $C = \partial R$ is the bounding curve, ${\bf r}(t) = x(t) \, {\bf i} + y(t) \, {\bf j}$ is a parametrization of $C$ in a counterclockwise direction and ${\bf F} = M \, {\bf i} + N {\bf j}$ is the field. To get the version using flux and divergence, first write the above (using simplified notation) as
$$\oint_{\partial R} \langle M, N \rangle \cdot \langle x'(t), y'(t) \rangle \, dt = \iint (N_x - M_y) \, dA$$.
Now consider the field ${\bf G} = \langle -N, M \rangle$. Applying Green's theorem to this field, we have that
$$\oint_{\partial R} \langle -N, M \rangle \cdot \langle x'(t), y'(t) \rangle, dt = \iint (M_x + N_y) \, dA.$$
The integrand of the line integral is also equal to $\langle M, N \rangle \cdot \langle y'(t), -x'(t) \rangle \, dt$, which is equal to ${\bf F} \cdot {\bf n} \, ds$ for the same reason that ${\bf F} \cdot {\bf T} \, ds = {\bf F} \cdot d{\bf r}$, where ${\bf T}$ is the unit tangent vector. Notice that ${\bf n}$ is obtained from ${\bf T}$ by rotating it counterclockwise by 90 degrees. Since the curve was traversed in a counterclockwise direction, the field ${\bf n}$ is outward pointing.
That the two integrals, $\oint_{\partial R} {\bf F} \cdot {\bf n} \, ds$ and $\iint_{R} {\rm div} \, {\bf F} \, dA$, are equal can be argued intuitively by saying that summing the divergence over interior points leads to a canceling out of inward/outward flow. What remains is the flow out of the boundary.
If you already have developed an intuition for the physics, then you can read Green's theorem as a 2D version of Gauss's law: the flux through the bounding curve is proportional to the total charge in the enclosed region.