We suggest a novel mathematical framework for the in-homogeneous spatial spreading of an infectious disease in human population, with particular attention to COVID-19. Common epidemiological models, e.g., the well-known susceptible-exposed-infectious-recovered (SEIR) model, implicitly assume uniform (random) encounters between the infectious and susceptible sub-populations, resulting in homogeneous spatial distributions. However, in human population, especially under different levels of mobility restrictions, this assumption is likely to fail. Splitting the geographic region under study into areal nodes, and assuming infection kinetics within nodes and between nearest-neighbor nodes, we arrive into a continuous, "reaction-diffusion", spatial model. To account for COVID-19, the model includes five different sub-populations, in which the infectious sub-population is split into pre-symptomatic and symptomatic. Our model accounts for the spreading evolution of infectious population domains from initial epicenters, leading to different regimes of sub-exponential (e.g., power-law) growth. Importantly, we also account for the variable geographic density of the population, that can strongly enhance or suppress infection spreading. For instance, we show how weakly infected regions surrounding a densely populated area can cause rapid migration of the infection towards the populated area. Predicted infection "heat-maps" show remarkable similarity to publicly available heat-maps, e.g., from South Carolina. We further demonstrate how localized lockdown/quarantine conditions can slow down the spreading of disease from epicenters. Application of our model in different countries can provide a useful predictive tool for the authorities, in particular, for planning strong lockdown measures in localized areas-such as those underway in a few countries.