TransWikia.com

Extract population count from GeoTIFF based on mask layer

Geographic Information Systems Asked on August 17, 2021

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

2 Answers

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)

Correct answer by snowman2 on August 17, 2021

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.

Answered by Jose on August 17, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP