Raster smoothing - assigning each raster the average value of its neighbors

Question

I have a raster world map of population densities at a resolution of 30 arc seconds. For each raster, I want to find out the average density of neighbors within a certain distance (e.g. 0.1 degrees). The aim is to keep and save in an extra file only those raster values with a density greater than a certain threshold (e.g. 1000 people per sqkm). All contiguous rasters above this threshold will then be clustered to form a core city. Next, I want to find the city fringe to those core cities. This shall be composed by rasters (attached to the core cities) with population densities between 500 and 1000. Core plus fringe will compose a complete city. Those cities then should be transformed to a polygon shape file where I can assign each city an own ID. Furthermore, having a polygon shape file allows me to easily assign values from other raster maps to the city using zonal statistics.

My question is how to do all of this. I tried around a little with Warp but I think that the result is different from what I want because Warp only is able to average over the neighbors to reproject the small to bigger raster, but which rasters are defined as neighbors depends on the layer extent chosen.

Below a sample image of my raster map.

Is there a way without using GRASS? Because I cannot get GRASS running..

Answer 1

If you're interested in using GRASS for this, here are the steps. To smooth the original, r.neighbors creates a new raster where each pixel gets the average value of surrounding pixels within a window. THen you can use r.recode to get a new raster with discrete categorical values replacing the continous population density. And finally, r.to.vect will convert to vector polygons. From the GRASS command shell:

g.region -ap raster=pop_dens
# Window size of 12 means 12 pixels of 30 arcsec = 6 mins = 0.1 degrees 
r.neighbors input=pop_dens output=pop_dens_smoothed method=average size=12
# Recode the continuous raster to only three values at the breaks you need    
r.recode input=pop_dens_smoothed output=pop_dens_recode rules=- <<EOF
0:500:0
501:1000:500
1000:99999:1000
EOF

r.to.vect input=pop_dens_recode output=pop_dens type=area column=density

Of course You can also run these commands from the GRASS GUI, or from the Processing Toolbox in QGIS.

I would mention that if your original population density is in a non-projected CRS (you mention resolution in degrees), then I wouldn't trust the densities at high latitudes.