Consider a sublattice $\langle(a_x,a_y),(b_x,b_y)\rangle$ of index $n$. A node $(r_x,r_y)$ belongs to that sublattice when $A\begin{pmatrix}k_a\\k_b\end{pmatrix}=\begin{pmatrix}r_x\\r_y\end{pmatrix}$ where $A=\begin{pmatrix}a_x&b_x\\a_y&b_y\end{pmatrix}$ and $k_{x,y}$ are some integers. Since $|\det A| = n$ (we can even assume $\det A=n$ without loss of generality), we have: $\begin{pmatrix}k_a\\k_b\end{pmatrix}=\displaystyle\frac{1}{n}\begin{pmatrix}b_y&-b_x\\-a_y&a_x\end{pmatrix}\begin{pmatrix}r_x\\r_y\end{pmatrix}$, so the solution of that linear equation system is surely integral if $r_{x,y}$ are both divisible by $n$.
In particular, the node $(n,0)$ is contained in every sublattice of index $n$. There may be other nodes between $(0,0)$ and $(n,0)$ if the distance $d$ between them is a divisor of $n$. Thus, for every $d|n$ there is a way to select one basis vector of the sublattice as $(d,0)$. Another basis vector can be chosen in the form $(a,n/d)$ so the index is indeed $d\times n/d-0\times a=n$, and two choices of $a$ lead to the same sublattice iff their difference is a multiple of $d$. Hence the formula $\sum_{d|n}d$ from your comment.
For the general $D$-dimensional case, it works by induction: for every divisor $d$ of $n$, we can select the first basis vector to be $(d,0,\ldots,0)$; and to select the remaining vectors, we look at every index $n/d$ sublattice in $(D-1)$-dimensional space and take it into account multiple times, namely the number of possible shifts of that $(D-1)$-dimensional sublattice that do not transform it to itself again, which equals to its index, $n/d$; it's the same number for all sublattices, so instead of summing over $(D-1)$-dimensional sublattices, we count them and multiply their number by $n/d$. After replacement $d\leftrightarrow n/d$ we get the recurrent formula from OEIS A128119/A160870: $T_{n,D} = \sum_{d|n}d\times T_{d,D-1}$ with $T_{1,D}=T_{n,1}=1$.