I am attempting to use the zonal() function of the Terra package in R to find out how much population there is for a given municipality. I am using nic_m_0_2020.tif as the source for the population data and OSM admin_level_6_nicaragua for the municipalities which i have verified has valid geometry and the correct CRS. I also tested using a different municipal boundary layer.


municipios <- vect("/home/benja/Downloads/Pop Nica/NIC.shp", )

onvres <- rast("/home/benja/Downloads/Pop Nica/nic_m_0_2020.tif")

mujeres <- rast("/home/benja/Downloads/Pop Nica/nic_f_0_2020.tif")

pop_onvres <- zonal(onvres, municipios, fun = sum, as.raster = T)

After I run this code, I end up with incorrect results.

When run in qgis, the result looks like this:

enter image description here

note that the max and minimum for my values calculated in R do not match the results from QGIS:

> print(pop_onvres)
class       : SpatRaster 
dimensions  : 5181, 6830, 1  (nrow, ncol, nlyr)
resolution  : 0.0008333333, 0.0008333333  (x, y)
extent      : -87.69125, -81.99958, 10.70792, 15.02542  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : spat_7HpJcOClpMquNmS_8284.tif 
name        : nic_m_0_2020 
min value   :     37.11481 
max value   :   1079.65405 

If I plot the population raster, I get this:

enter image description here

When I overlay my zonal statitics layer, I get this:

enter image description here

You can see that there is not a generated zonal polygon for every polygon in my shapefile, and there are places in the population raster that have data points that were not included in the zonals.

The full country municipal layer looks like this when overlayed on the population raster.:

enter image description here

My result only includes this subset of the polygons:

enter image description here

I cannot seem to sipmly get the zonal() function to just calculate the sum of population per polygon. It should be a fairly straightforward and basic operation and I cannot get it to work.

It does seem to be counting correctly, but only doing it for a handful of my polygons.

  • 3
    You likely have NA's in your raster. The sum function will return an NA if there are any NA's in the vector. Try adding the argument na.rm=TRUE after fun=sum and see if that makes any difference. Commented Nov 6, 2023 at 22:31
  • That seems to have done did it. Commented Nov 7, 2023 at 0:19

1 Answer 1


As per comment by Jeffrey Evans, my raster has NAs and adding na.rm = T to the zonal function fixed it.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.