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 · 2020-02-01T01:07:14.293

3

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)

edited Feb 01 '20 at 01:07

answered Jan 27 '20 at 22:18

snowman2

7,321
12
29
54

Thank you for your answer! Unfortunately, I was not able to reproduce it yet due to an error with GDAL CRSError: The EPSG code is unknown.. The error is described here https://rasterio.readthedocs.io/en/latest/faq.html but I was not able to solve it yet. – Stücke Jan 29 '20 at 07:25
Do you have the environment variable GDAL_DATA set? Also is PROJ_LIB set? – snowman2 Jan 29 '20 at 13:52
To be honest, I don't even know what that means. – Stücke Jan 29 '20 at 14:13
What operating system are you using? – snowman2 Jan 29 '20 at 15:05
I am using Windows 10. – Stücke Jan 30 '20 at 15:46
How did you install it? – snowman2 Jan 30 '20 at 17:42
I believe I just installed it either via pip install gdal or conda install -c conda-forge gdal. – Stücke Jan 31 '20 at 12:22
Just updated installation method. It should work with conda-forge. If you activate your environment before use, it should set the environment variables for you. – snowman2 Feb 01 '20 at 01:08
Thank you for cour continuous support. I deinstalled gdal and reinstalled it via conda install -c conda-forge gdal. I was able to successfully install it without any error but I do receive a new error when I run the script https://pastebin.com/vDSXxmGk – Stücke Feb 03 '20 at 07:24
Looks like UBA isn't in your environment variables. Not sure what that is. – snowman2 Feb 04 '20 at 02:29
It's the name of my environment. Strange ... – Stücke Feb 04 '20 at 07:22

score 1 · Answer 2 · answered Apr 27 '20 at 15:35

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 Answers2