Optimize station coverage

I have ~10000 irregularly placed nodes in x,y space. Each node must be covered by a station. A node is defined as covered by a station if it is within 1000m of the station. Each station can only cover up to 4 nodes, regardless of how many nodes fall within 1000m of the station.

I want to optimize the location of the stations such that the total number of stations is minimized with the constraint that every node must be covered. The algorithm needs to be able to manipulate the x,y location of the stations, as well as add and remove stations. The add and remove stations part is the piece that I am not sure how to implement.

Looking for advice on terms/algorithms that can help me get started. I'm guessing similar optimizations are conducted for placing cell and wifi towers?, but I can find anything. I intend to implement the solution in python with scikit-learn