0

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):

enter image description here

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

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,
)

enter image description here

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)
11
  • 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 rasterio.readthedocs.io/en/latest/faq.html but I was not able to solve it yet.
    – Stücke
    Jan 29 '20 at 7: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
1

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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