Percent of polygon population / area experiencing deficit/surplus

example.sh

#!/usr/bin/env bash
set -e

YEARMON=$1
# Computes the percent of land area and population of each country that is experiencing water anomalies (surplus or deficit) 
# of return periods >= 10 or 30 years. Output as two CSVs (one for population percent and one for area percent) each with
# four coulums: (pct exposed to 10-yr deficit, pct exposed to 30-yr deficit, pct exposed to 10-yr surplus, pct exposed to 30-yr surplus)
# Uses exactextract to handle pixels that are partially covered by country polygons.
# Uses gdal_calc to make a set of four mask rasters that are passed as weights (0 or 1) to exactextract.

#GPW="/home/dbaston/Downloads/gpw-v4-population-count-rev10_2015_30_sec_tif/gpw_v4_population_count_rev10_2015_30_sec.tif"
#LAND="/home/dbaston/Downloads/gpw-v4-land-water-area-rev10_landareakm_30_sec_tif (3)/gpw_v4_land_water_area_rev10_landareakm_30_sec.tif"
COUNTRIES="/home/dbaston/Downloads/ne_10m_admin_0_countries.shp"

# I couldn't use the GPW population count or land area rasters directly because the origin point of the grid is incompatible with the
# surplus and deficit data I was using.
# gdalinfo reports the GPW origin point as:
# Origin = (-180.000000000000000,84.999999999991758)
# In the "fixed" versions in /tmp/ below, I fudged the origin point over to (-180, 85).
GPW="/tmp/gpw_fixed.tif"
LAND="/tmp/landarea_fixed.tif"

WSIM_COMPOSITE_DIR="/mnt/fig/WSIM/WSIM_derived_V2/composite_adjusted"

DEFICIT_10=`tempfile -s ".tif"`
DEFICIT_30=`tempfile -s ".tif"`
SURPLUS_10=`tempfile -s ".tif"`
SURPLUS_30=`tempfile -s ".tif"`
gdal_calc.py -A "NETCDF:${WSIM_COMPOSITE_DIR}/composite_adjusted_12mo_${YEARMON}.nc:deficit" --calc "A<=-10" --type=Byte --outfile ${DEFICIT_10} --overwrite --co "COMPRESS=DEFLATE"
gdal_calc.py -A "NETCDF:${WSIM_COMPOSITE_DIR}/composite_adjusted_12mo_${YEARMON}.nc:deficit" --calc "A<=-30" --type=Byte --outfile ${DEFICIT_30} --overwrite --co "COMPRESS=DEFLATE"
gdal_calc.py -A "NETCDF:${WSIM_COMPOSITE_DIR}/composite_adjusted_1mo_${YEARMON}.nc:surplus"  --calc "A>=10"  --type=Byte --outfile ${SURPLUS_10} --overwrite --co "COMPRESS=DEFLATE"
gdal_calc.py -A "NETCDF:${WSIM_COMPOSITE_DIR}/composite_adjusted_1mo_${YEARMON}.nc:surplus"  --calc "A>=30"  --type=Byte --outfile ${SURPLUS_30} --overwrite --co "COMPRESS=DEFLATE"

OUTFILE_POP="/tmp/pct_pop_experiencing_deficit_${YEARMON}.csv"
OUTFILE_LAND="/tmp/pct_land_experiencing_deficit_${YEARMON}.csv"

exactextract -p ${COUNTRIES} -f NAME -r "${GPW}" -w ${DEFICIT_10} -w ${DEFICIT_30} -w ${SURPLUS_10} -w ${SURPLUS_30} -s "weighted fraction" -o "${OUTFILE_POP}" --progress
exactextract -p ${COUNTRIES} -f NAME -r "${LAND}" -w ${DEFICIT_10} -w ${DEFICIT_30} -w ${SURPLUS_10} -w ${SURPLUS_30} -s "weighted fraction" -o "${OUTFILE_LAND}"  --progress