Preliminary update to water demand historic tool dataset instructions…

… and scripts

requirements.txt

@@ -31,6 +31,7 @@ python-dateutil>=2.8.2
 pytz>=2023.3
 PyYAML>=6.0
 rasterio>=1.2.10
+scikit-image>=0.19.3
 scipy>=1.7.3
 six>=1.16.0
 tomli>=2.0.1

src/lisfloodutilities/water-demand-historic/README.md

@@ -9,18 +9,19 @@ The global sectoral water demand maps at 3 arcmin (or 0.05 degrees) resolution,
 # Data
 
 The utility consists of five scripts which require several following datasets and files. The full list of datasets is provided below. For each dataset, this documentation provides the link and the instructions relative to available version at the time of creation of the scripts. The correct functioning of the scripts is guaranteed for those datasets. Clearly, updated versions are regularly published: one relevant example is [Global Human Settlement - Datasets - European Commission (europa.eu)](https://ghsl.jrc.ec.europa.eu/datasets.php) (see Global Human Settlement Layer (GHSL)). Albeit the functioning of the utility is expected to be maintained, the users are invited to report problems or suggestions via GitHub issues.
-1. [FAO AQUASTAT](http://fao.org/aquastat/statistics/query/index.html?lang=en) sectoral water withdrawal estimates. Select "All Countries" and the following seven fields in the variable groups "Geography and population" and "Water use": "Gross Domestic Product (GDP)", "Industry, value added to GDP", "Agricultural water withdrawal", "Industrial water withdrawal", "Municipal water withdrawal", "Total water withdrawal", and "Irrigation water withdrawal". Select All Years, then click on the "Show Data" button.  On the new page, click on the "Download" button that will appear on the top right corner to download the xlsx file containing all data. Convert the Excel file to CVS and save it as `aquastat_clean.csv` in `aquastat_folder` specified in the configuration file.
-1. [Global Human Settlement Layer (GHSL)](https://ghsl.jrc.ec.europa.eu/ghs_pop2019.php) POP R2019A residential population estimates for target years 1975, 1990, 2000, and 2015. You will find data as an archive [here]( https://jeodpp.jrc.ec.europa.eu/ftp/jrc-opendata/GHSL/GHS_POP_MT_GLOBE_R2019A/). For each year, download the 0.0025° (9 arcsec, folders ending by "_9ss", e.g. "GHS_POP_E2015_GLOBE_R2019A_4326_9ss") maps in WGS84 projection and resample them to 0.01° using `gdalwarp -t_srs EPSG:4326 -tr 0.01 0.01 -r average -of GTiff GHS_POP_E<year>_GLOBE_R2019A_4326_9ss_V1_0.tif GHS_POP_E<year>_GLOBE_R2019A_4326_9ss_V1_0_reprojected.tif`, where `<year>` is 1975, 1990, 2000, or 2015. We resample the maps using `average` instead of `sum` as the latter is broken. Put resampled files in `ghsl_folder`.
-1. [Global Change Analysis Model (GCAM)](https://github.com/JGCRI/gcam-core/releases) regional water withdrawal and electricity consumption estimates. Download the Windows release package version 5.4, execute `run-gcam.bat`, wait for the model to finish. Execute `run-model-interface.bat`, click "File" > "Open" > "DB Open", select `output/database_basexdb`, and select all scenarios and all regions. Select `water demand`, select `water withdrawals by sector`, and click "Run query". Select `energy transformation`, select `electricity`, select `elec consumption by demand sector`, and click "Run query". For each tab, select all data with ctrl-a and click "Edit" > "Copy". Open a blank [Google Sheets](http://sheets.google.com/) spreadsheet, press ctrl-v, manually add the headers, click "File" > "Download" > "Comma-separated values", save as `water_demand.csv` and `elec_consumption.csv`, respectively, and put both files in `gcam_folder`.
+1. [FAO AQUASTAT](http://fao.org/aquastat/statistics/query/index.html?lang=en) sectoral water withdrawal estimates. Select "All Countries". Then Select the variable groups "Geography and population" and "Water use" (if needed select the subgroubps "Economy, development and food security" and "Water withdrawal by sector"). From these variable groups, select the following seven fields: "Gross Domestic Product (GDP)", "Industry, value added to GDP", "Agricultural water withdrawal", "Industrial water withdrawal", "Municipal water withdrawal", "Total water withdrawal", and "Irrigation water withdrawal". Select All Years, or manually select the years from 1976 to 2021, then click on the "Show Data" button if needed. Finally click on the "Download" button on the top right corner to download the xlsx or csv file containing all data. Convert the Excel file to CVS if needed and save it as `aquastat_clean.csv` in `aquastat_folder` specified in the configuration file.
+1. [Global Human Settlement Layer (GHSL)](https://human-settlement.emergency.copernicus.eu/ghs_pop2023.php) POP R2023A residential population estimates for target years from 1975 to 2020. You will find data as an archive [here](https://jeodpp.jrc.ec.europa.eu/ftp/jrc-opendata/GHSL/GHS_POP_GLOBE_R2023A/). For each year, download the 3 arcsec (folders ending by "_3ss"), e.g. "GHS_POP_E2015_GLOBE_R2023A_4326_3ss") maps in WGS84 projection and resample them to 0.01° using `gdalwarp -t_srs EPSG:4326 -tr 0.01 0.01 -r average -of GTiff GHS_POP_E<year>_GLOBE_R2023A_4326_3ss_V1_0.tif GHS_POP_E<year>_GLOBE_R2023A_4326_3ss_V1_0_reprojected.tif`, where `<year>` is 1975 to 2020. We resample the maps using `average` instead of `sum` as the latter is broken. Put resampled files in `ghsl_folder`. You will need to specify also the 'ghsl_folder_factor = 144' to account for resampling from 0.000833 (3 arcsec) to 0.01 degrees.
+1. [Global Change Analysis Model (GCAM)](https://github.com/JGCRI/gcam-core/releases) regional water withdrawal and electricity consumption estimates. Download the Windows release package version 7.1, execute `run-gcam.bat`, wait for the model to finish. Execute `run-model-interface.bat`, click "File" > "Open" > "DB Open", select `output/database_basexdb`, and select all scenarios and all regions. Select `water demand`, select `water withdrawals by sector`, and click "Run query". Select `energy transformation`, select `electricity`, select `elec consumption by demand sector`, and click "Run query". For each tab, select all data with ctrl-a and click "Edit" > "Copy". Open a blank [Google Sheets](http://sheets.google.com/) spreadsheet, press ctrl-v, manually add the headers, click "File" > "Download" > "Comma-separated values", save as `water_demand.csv` and `elec_consumption.csv`, respectively, and put both files in `gcam_folder`.
 1. [Gridded Livestock of the World (GLW 3)](https://doi.org/10.1038/sdata.2018.227) species distribution dataset. Download the eight zip files (each representing a different species) and extract them to `glw_folder`.
 1. [Huang et al. (2018)](https://doi.org/10.5281/zenodo.1209296) global gridded water withdrawal estimates. These estimates are not incorporated in our dataset, but are only used for the sake of comparison. Download all 7z files and extract them to `huang_folder`.
 1. [Multi-Source Weather (MSWX)](http://www.gloh2o.org/mswx) daily and monthly mean air temperature. Download using rclone as explained in the FAQ on the web page. Put the daily and monthly netCDFs in `<mswx_folder>/Past/Temp/Daily` and `<mswx_folder>/Past/Temp/Monthly`, respectively, where `mswx_folder` is specified in the configuration file.
-1. [US Census Bureau](https://www.census.gov/geographies/mapping-files/time-series/geo/carto-boundary-file.html) state borders shape file. Download `cb_2018_us_state_500k.zip`, unzip it, open the shape file in QGIS, open the Field Calculator, enter `STATEFP` in "Expression", enter `STATEFP_` in "Output field name", and click OK. Rasterize to 0.01° using `gdal_rasterize -l cb_2018_us_state_500k -a STATEFP_ -tr 0.01 0.01 -a_nodata 0.0 -te -180.0 -90.0 180.0 90.0 -ot Float32 -of GTiff cb_2018_us_state_500k.shp cb_2018_us_state_500k_rasterized.tif`. Put `cb_2018_us_state_500k_rasterized.tif` in `us_states_folder`.
+1. [US Census Bureau](https://www.census.gov/geographies/mapping-files/time-series/geo/cartographic-boundary.html) state borders shape file. Download `cb_2023_us_state_500k.zip`, unzip it, open the shape file in QGIS, open the Field Calculator, enter `STATEFP` in "Expression", enter `STATEFP_` in "Output field name", and click OK. Rasterize to 0.01° using `gdal_rasterize -l cb_2023_us_state_500k -a STATEFP_ -tr 0.01 0.01 -a_nodata 0.0 -te -180.0 -90.0 180.0 90.0 -ot Float32 -of GTiff cb_2023_us_state_500k.shp cb_2023_us_state_500k_rasterized.tif`. Put `cb_2023_us_state_500k_rasterized.tif` in `us_states_folder`.
 1. [USGS NWIS](https://waterdata.usgs.gov/nv/nwis/wu) water withdrawal estimates for 1985--present. For each state in the "Geographic Area" drop-down menu, select "State Data", "ALL Years", "State Total", and "ALL Categories" and click "Submit". Then select "Tab-separated data" and click "Submit". Do this for each state and put files in `usgs_water_use_folder`.
 1. USGS water withdrawal data files for [1985](https://water.usgs.gov/watuse/data/1985/index.html) and [1990](https://water.usgs.gov/watuse/data/1990/index.html) (to supplement the NWIS data). Download "Data file of state-level data" for each year (`us85st.txt` and `us90st.txt`) and put the files in `usgs_water_use_folder`.
 1. [Vassolo and Döll (2005)](https://doi.org/10.1029/2004WR003360) industrial and thermoelectric water withdrawal maps (included in this repository with permission from Petra Döll). The industrial withdrawal shape file is rasterized using `gdal_rasterize -l "industry_paper_vassolo&doell" -a MANUF_WWD -tr 0.5 0.5 -a_nodata 0.0 -te -180.0 -90.0 180.0 90.0 -ot Float32 -of GTiff "industry_paper_vassolo&doell.shp" manuf_wwd.tif`. Thermoelectric shape file rasterized using same command but with `WWD_PS` and `wwd_ps.tif`. Put the files in `vassolo_doll_folder`.
-1. [Thematic Mapping](https://thematicmapping.org/downloads/world_borders.php) country borders shape file. Download `TM_WORLD_BORDERS-0.3.zip`, unzip, and rasterize to 0.01° using `gdal_rasterize -l TM_WORLD_BORDERS-0.3 -a UN -tr 0.01 0.01 -a_nodata 0.0 -te -180.0 -90.0 180.0 90.0 -ot Float32 -of GTiff TM_WORLD_BORDERS-0.3.shp TM_WORLD_BORDERS_UN_rasterized.tif`. Put the result in `world_borders_folder`.
-1. [World Bank](https://data.worldbank.org/) manufacturing value added and gross domestic product data. Search for "Manufacturing, value added (constant 2010 US$)" and "GDP (constant 2010 US$)", download as CSV, and put in `world_bank_folder` (Remove Metadata csv files, if any).
+1. [GISCO](https://ec.europa.eu/eurostat/web/gisco/geodata/administrative-units/countries) country borders shape file. Select year 2024, File format SHP, Geometry type Polygons (RG), Scale 01M, Coordinate system EPGS 4326, click Download. Rasterize to 0.01° using `gdal_rasterize -l CNTR_RG_01M_2024_4326 -a UN -tr 0.01 0.01 -a_nodata 0.0 -te -180.0 -90.0 180.0 90.0 -ot Float32 -of GTiff CNTR_RG_01M_2024_4326.shp CNTR_RG_01M_2024_4326_rasterized.tif`. Put the result in `world_borders_folder`.
+
+1. [World Bank](https://data.worldbank.org/) manufacturing value added and gross domestic product data. Search for "Manufacturing, value added (constant 2015 US$)" and "GDP (constant 2015 US$)", download as CSV, and put in `world_bank_folder` (Remove Metadata csv files, if any).
 <!---
 1. USGS Circular xxx 1980 (https://pubs.usgs.gov/circ/1983/1001/report.pdf) Table 7 provided in ancillary_data folder
 1. [Global Power Plant Database](https://datasets.wri.org/dataset/globalpowerplantdatabase)

src/lisfloodutilities/water-demand-historic/config_global_3arcmin.cfg

@@ -0,0 +1,23 @@
+year_start = 1975
+year_end = 2023
+templatemap_path = /media/sf_VMSharedFolder/WaterDemandMapsDataset2024/template_Global_03min.nc
+output_folder = /media/sf_VMSharedFolder/WaterDemandMapsDataset2024/vmtest_lisflood-water-demand-historic_global_3arcmin
+ancillary_data_folder = ancillary_data
+aquastat_folder = /media/sf_VMSharedFolder/WaterDemandMapsDataset2024/AQUASTAT
+#eia_folder = /eos/jeodpp/data/projects/WEFE/shared/EIA_electricity_database
+gcam_folder = /media/sf_VMSharedFolder/WaterDemandMapsDataset2024/GCAM_v7.1
+ghsl_folder = /media/sf_VMSharedFolder/WaterDemandMapsDataset2024/GHS
+ghsl_folder_factor = 144
+glw_folder = /media/sf_VMSharedFolder/WaterDemandMapsDataset2024/GLW
+huang_folder = /media/sf_VMSharedFolder/WaterDemandMapsDataset2024/HUANG
+#htap_folder = /eos/jeodpp/data/projects/WEFE/shared/OMI-HTAP
+#lohrmann_power_plant_db_folder = /eos/jeodpp/data/projects/WEFE/shared/lohrmann_power_plant_database
+mswx_folder = /media/sf_VMSharedFolder/WaterDemandMapsDataset2024/MSWX
+#pku_fuel_folder = /eos/jeodpp/data/projects/WEFE/shared/PKU-FUEL
+#un_icsd_folder = /eos/jeodpp/data/projects/WEFE/shared/UN_Industrial_Commodity_Stats
+us_states_folder = /media/sf_VMSharedFolder/WaterDemandMapsDataset2024/US_STATES
+usgs_water_use_folder = /media/sf_VMSharedFolder/WaterDemandMapsDataset2024/USGS
+vassolo_doll_folder = industrial_water_use_data_Vassolo_Doll
+world_borders_folder = /media/sf_VMSharedFolder/WaterDemandMapsDataset2024/WORLDBORDERS
+world_bank_folder = /media/sf_VMSharedFolder/WaterDemandMapsDataset2024/WORLDBANK
+#wri_power_plant_db_folder = /eos/jeodpp/data/projects/WEFE/shared/WRI_power_plant_database

src/lisfloodutilities/water-demand-historic/step1_population_density.py

@@ -18,7 +18,7 @@
 import pandas as pd
 import numpy as np
 from netCDF4 import Dataset
-import matplotlib.pyplot as plt
+#import matplotlib.pyplot as plt
 from tools import *
 import rasterio
 
@@ -74,9 +74,9 @@ def main():
         print('Loading and upscaling '+str(pop_years[ii])+' population data')
 
         # Load raw data (total population per grid-cell)
-        # Factor 16 to account for resampling from 0.0025 to 0.01 degrees
-        src = rasterio.open(os.path.join(config['ghsl_folder'],'GHS_POP_E'+str(pop_years[ii])+'_GLOBE_R2019A_4326_9ss_V1_0','GHS_POP_E'+str(pop_years[ii])+'_GLOBE_R2019A_4326_9ss_V1_0_reprojected.tif'))
-        pop_data = src.read(1).astype(np.single)*16
+        # Factor ghsl_folder_factor = 144 to account for resampling from 0.000833 (3 arcsec) to 0.01 degrees
+        src = rasterio.open(os.path.join(config['ghsl_folder'],'GHS_POP_E'+str(pop_years[ii])+'_GLOBE_R2023A_4326_3ss_V1_0','GHS_POP_E'+str(pop_years[ii])+'_GLOBE_R2023A_4326_3ss_V1_0_reprojected.tif'))
+        pop_data = src.read(1).astype(np.single)*config['ghsl_folder_factor']
         pop_data[pop_data<0] = 0
         refmatrix = src.get_transform()
         src.close()

src/lisfloodutilities/water-demand-historic/step2_domestic_demand.py

@@ -98,7 +98,7 @@ def main():
     print('-------------------------------------------------------------------------------')
     print('Loading and resampling country border raster')
     t0 = time.time()
-    country_code_map = load_country_code_map(os.path.join(config['world_borders_folder'],'TM_WORLD_BORDERS_UN_rasterized.tif'),mapsize_global)
+    country_code_map = load_country_code_map(os.path.join(config['world_borders_folder'],'CNTR_RG_01M_2024_4326_rasterized.tif'),mapsize_global)
     print("Time elapsed is "+str(time.time()-t0)+" sec")
 
 

src/lisfloodutilities/water-demand-historic/step3a_industrial_demand.py

@@ -71,7 +71,7 @@ def main():
     print('-------------------------------------------------------------------------------')
     print('Loading and resampling country border raster')
     t0 = time.time()
-    country_code_map = load_country_code_map(os.path.join(config['world_borders_folder'],'TM_WORLD_BORDERS_UN_rasterized.tif'),mapsize_global)
+    country_code_map = load_country_code_map(os.path.join(config['world_borders_folder'],'CNTR_RG_01M_2024_4326_rasterized.tif'),mapsize_global)
     country_code_map_360x720 = resize(country_code_map,(360,720),order=0,mode='edge',anti_aliasing=False)
     print("Time elapsed is "+str(time.time()-t0)+" sec")
 

src/lisfloodutilities/water-demand-historic/step3b_industrial_demand.py

@@ -71,7 +71,7 @@ def main():
     print('-------------------------------------------------------------------------------')
     print('Loading and resampling country border raster')
     t0 = time.time()
-    country_code_map = load_country_code_map(os.path.join(config['world_borders_folder'],'TM_WORLD_BORDERS_UN_rasterized.tif'),mapsize_global)
+    country_code_map = load_country_code_map(os.path.join(config['world_borders_folder'],'CNTR_RG_01M_2024_4326_rasterized.tif'),mapsize_global)
     country_code_map = fill(country_code_map)
     country_code_map_1800x3600 = resize(country_code_map,(1800,3600),order=0,mode='edge',anti_aliasing=False)
     print("Time elapsed is "+str(time.time()-t0)+" sec")

src/lisfloodutilities/water-demand-historic/step4_livestock_demand.py

@@ -77,7 +77,7 @@ def main():
     print('-------------------------------------------------------------------------------')
     print('Loading and resampling country border raster')
     t0 = time.time()
-    country_code_map = load_country_code_map(os.path.join(config['world_borders_folder'],'TM_WORLD_BORDERS_UN_rasterized.tif'),mapsize_global)
+    country_code_map = load_country_code_map(os.path.join(config['world_borders_folder'],'CNTR_RG_01M_2024_4326_rasterized.tif'),mapsize_global)
     country_code_map = fill(country_code_map)
     country_code_map_1800x3600 = resize(country_code_map,(1800,3600),order=0,mode='edge',anti_aliasing=False)
     print("Time elapsed is "+str(time.time()-t0)+" sec")