Set Cover
If you restrict the possible center points of the circles to be the set of points, then this is an instance of set cover where the universe has $|P| \le 1000$ elements and where each set contains $N$ elements. While set cover is NP-hard, there are algorithms for set cover that are faster than exhaustive search.
For these parameters, I expect it should be possible to find an optimal or near-optimal solution to the associated set cover problem, using standard algorithms for set cover. I would recommend starting with the integer linear programming formulation and using a standard ILP solver (Gurobi, CPLEX) to find an optimal or near-optimal solution to the set cover problem. That might be efficient enough for your purposes.
Your approach is known as the greedy algorithm for set cover. It is known to give a $\lg N$-factor approximation, i.e., a covering that uses no more than $\lg N$ times as many circles as the optimal covering. Thus, the ILP formulation may yield better solutions. There are also heuristics for improving the greedy algorithm, e.g., BubbleSearch, which makes small random tweaks to the order suggested by the greedy strategy, repeats that some number of times, and takes the best solution found so far. You could try them and compare them to how well the ILP formulation works, if you wish.
Geometric Set Cover
You could also look at algorithms for geometric set cover; some of those techniques may be useful here. In particular, there are approximation algorithms and fast heuristics for geometric set cover, and you could try seeing if they are useful in your particular setting. They use the extra geometric structure provided by the problem: this is not a totally arbitrary instance of set cover, but there is some additional structure.
Arbitrary Center Points
The above algorithms have one critical limitation: they only consider circles that are centered at some point in $P$. There's no reason to expect the optimal solution to your problem to have this property. Thus, they might fail to find the optimal solution to your original problem.
Fortunately, it's possible to fix this limitation. It's possible to come up with a not-too-large set $U$ of candidate circles, such that the optimal solution is obtained by taking a subset of circles from $U$; then we can apply set-cover techniques to pick a minimal-size subset of $U$ that covers all of $P$. In particular, we will use the following geometric facts:
Lemma 1. Let $C_0$ be any circle that covers exactly $N$ points. Let $S \subset P$ be the set of $N$ points covered by $C_0$. Then there exists a circle $C$ that (i) covers exactly $S$ and no other points and either (ii.1) has at least three points of $S$ along its circumference or (ii.2) has two points of $S$ along its circumference and is centered at their midpoint.
Proof. Consider any circle that covers exactly $S$ (and nothing more). Pick the point $p \in S$ that is closest to the boundary of this circle. Reduce the radius of the circle until $p$ is on the boundary. Now if this circle doesn't have two points of $S$ along its circumference yet, move the center of the circle towards $p$, reducing the radius at the same rate, until the first time some other point is on the boundary -- call it $q$. At this point we have a circle $C$ that covers all of $S$, contains $p,q$ on its boundary, and doesn't cover any points other than $S$. Now let $m$ be the midpoint between $p,q$, and let $L$ be the line that is perpendicular to $pq$ and that intersects $m$. If $C$ is centered at $m$, we're done. Otherwise, $C$'s center is somewhere on $L$. Move the center of $C$ along $L$ towards $m$, decreasing its radius at the same time so that $p,q$ remain on its circumference, until either the center reaches $m$ or the circumference hits a third point. Either way, we are done. $\Box$
Lemma 2. Given any three non-colinear points $p,q,r$, there is a unique circle that is tangent to all three (i.e., has all three on its circumference). This circle can be efficiently computed.
Proof. Use the same procedure described in the proof of Lemma 1, but now without worrying about how many points are covered by the circle. $\Box$
Lemma 3. Given any two points $p,q$, there is a unique circle centered at their midpoint and tangent to both $p,q$.
Proof. Obvious. $\Box$
With these lemmas, we can now suggest an algorithm:
For each set of three points $p,q,r \in P$, find the circle that is tangent to $p,q,r$. Identify the points covered by this circle (e.g., using a quadtree). If this circle covers exactly $N$ points, add it to $U$.
For each set of two points $p,q \in P$, find the circle centered at their midpoint and tangent to both. If this circle covers exactly $N$ points, add it to $U$.
Lemma 1 now ensures that we need only consider circles in $U$; any set of $N$ points that can be exactly covered by some circle, can be covered by some circle in $U$. Therefore, there exists an optimal solution to the original problem that uses only circles in $U$. As an optimization, if there are two circles in $U$ that cover the same set of points, you can keep one of them and discard the other.
Therefore, the problem becomes:
Given a set $P$ of points and a set $U$ of circles (each covering exactly $N$ points from $P$), find the smallest subset of $U$ that covers all of $P$.
This is a finite set cover problem, and the solution to this set-cover problem is guaranteed to be optimal. We obtain a set-cover instance with a universe of size $|P| \le 1000$ and with $|U|$ sets, each of size $N$. This set cover problem can then be solved using standard algorithms for set cover, e.g., using an ILP solver. We'll have $|U| \le |P|^3/6$ at worst, though in practice I expect $|U|$ will typically be much smaller (closer to $|P|^2$ or so). The running time to create the set-cover instance will be at most $O(|P|^3)$; then we have to solve the set-cover instance. I anticipate that, for these parameters, it might be possible to obtain an optimal or very-nearly-optimal solution to the original problem using off-the-shelf ILP solvers.
My thanks to @j_random_hacker for substantial assistance with this answer.
