Extract population count from GeoTIFF based on mask layer

Question

I downloaded Population Count data from Gridded Population of the World (GPW), v4 as a *.GeoTIFF file. In a next step, I would like to extract the population count of a given area based on a mask layer (e.g. the political boundaries of Portugal imported from a *.geojson).

I came across an answer on StackExchange, which I believe is trying to achieve a similar result. However, I could not reproduce it due to NotImplementedError: multi-dimensional sub-views are not implemented and I suspect the solution would be overengineered for me.

When I open the *.GeoTIFF file in QGIS and display the layer properties, I see the following information (I suppose, I am interested in the STATISTICS_MEAN):

The files are avialble for download here (only for seven days or 100 downloads though):

GeoTIFF: https://send.firefox.com/download/a4be075f0dda3acb/#UlUv5RDx8i3XXvV4dQHtEA
Geojson: https://send.firefox.com/download/54d3458a17ca355a/#DM8M3hGOfcoIDOVeKcMZJA

snowman2 · Accepted Answer

This is something you can do with geocube and rioxarray.

Here is an example I followed to generate this:
https://corteva.github.io/geocube/stable/examples/zonal_statistics.html

import geopandas
import numpy
import rioxarray
from geocube.api.core import make_geocube

gpd = geopandas.read_file("portugal.geojson")
gpd["mask"] = 1

raster = rioxarray.open_rasterio(
    "gpw_v4_population_count_rev11_2020_30_sec_portugal.tif",
    mask_and_scale=True,
)
mask = make_geocube(
    gpd,
    measurements=["mask"],
    like=raster,
    fill=0,
)

Then, to get the total population:

total_population = raster.where(mask.mask).sum().item()

16079493.21870717

EDIT:

Installation method:

conda create -n geo geocube
conda activate geo
(geo)

Jose · Answer

You can do that directly with GDAL, which is likely installed if you have QGIS already, both in Python and in the console. Python version here:

import numpy as np
import matplotlib.pyplot as plt
import gdal

fname = "gpw_v4_population_count_rev11_2020_30_sec_portugal"
world = "TM_WORLD_BORDERS_SIMPL-0.3.shp"
data = gdal.Warp("", fname, format="MEM", 
                 cutlineDSName=world, cropToCutline=True,
                 dstNodata=-999, 
                 cutlineWhere="NAME='Portugal'"
                ).ReadAsArray()
print(f"Total population {data[data != -999].sum()}")

A couple of things:
1. I have cropped to only consider the envelope of Portugal
2. I have set the no data value to -999
3. Your links aren't working anymore, so I just 
    * Downloaded the original dataset (see filename)
    * Used the TM_WORLD_BORDERS_SIMPL dataset
4. From the TM_WORLD_BORDERS_SIMPL, I just select the Portugal entry.

Extract population count from GeoTIFF based on mask layer

2 Answers

Add your own answers!

Ask a Question