If city-dwelling versus non-city-dwelling is dichotomous, I see no problem with using that distinction as a categorical predictor. If you think that specific Provinces might differ in terms of baseline measures of your outcomes, or in terms of the relation of city/non-city-dwelling to differences in outcomes, you could include the Province as a random effect in your model. The Province for each individual could be included as a random effect for either or both of the intercept (baseline outcome values, say, for non-city) or the slope (difference between city and non-city-dwelling outcome values) in the model.
It's not clear, however, how you are dealing with locations like suburbs that contain some aspects of both city and non-city; I live in such a suburban town immediately next to a large city. Instead of a strict dichotomy, I wonder if you would be better off using something like a measure of local population density as a continuous predictor instead.
I'm assuming here that Provinces are somewhat equivalent to States in the US, and that each city dweller is also a resident of one particular Province.