Merge pull request GEMScienceTools#6 from GEMScienceTools/no_raster_i…

…n_boundary

Handle no raster in boundary cases + extend download_worldpop

.gitignore

@@ -1,3 +1,6 @@
+# Avoid large files
+data/worldpop/*.tif
+
 # Byte-compiled / optimized / DLL files
 __pycache__/
 *.py[cod]

README.md

@@ -16,21 +16,24 @@ The dependencies are as follows:
 
 * numpy
 * pandas
+* fiona
 * geopandas
 * rasterio
 * shapely^
 * osgeo
 * GDAL^^
 
-*^There is currently a known issue where importing from shapely returned an AssertionError when loading the GEOS library. This can be resolved by installing shapely before fiona, rasterio, and GDAL. See [this link](https://sgillies.net/2019/06/23/fix-for-geos-dll-bug-shapely-1-7a2.html) for more details. If that doesn't work, try using the command ``pip install shapely --no-binary shapely``*
+*^There is currently a known issue where importing from shapely returned an AssertionError when loading the GEOS library. This can be resolved by installing shapely before fiona, rasterio, and GDAL. See [this link](https://sgillies.net/2019/06/23/fix-for-geos-dll-bug-shapely-1-7a2.html) or [this link](https://github.com/Toblerity/Shapely/issues/553) for more details. If that doesn't work, try using the command ``pip install shapely --no-binary shapely``*
 
 
-*^^GDAL requires installation prior to ``pip`` installation. This can be done using ``brew``. Windows users might consider installing GDAL using [OSGeo4W](https://trac.osgeo.org/osgeo4w/). macOS users might consider using the [KyngChaos installer](https://www.kyngchaos.com/software/frameworks/)*
+*^^GDAL requires installation prior to ``pip`` installation. This can be done using ``brew``. Windows users might consider installing GDAL using [OSGeo4W](https://trac.osgeo.org/osgeo4w/). macOS users might consider using the [KyngChaos installer](https://www.kyngchaos.com/software/frameworks/). Additionally, if the ``pip`` installation fails, be sure to check that the versions between ``brew`` and ``pip`` correspond to one another.*
 
 ## Getting started
 
 First, ensure you are in the appropriate working directory (i.e. the main directory)
 
+### Downloading WorldPop data
+
 We will need raster data from WorldPop to successfully execute the demo. This file was not included with the repository due to the large file size.
 
 WorldPop data can be downloaded using the function **download_worldpop.py**.. This downloads the dataset to the expected location by **main_script.py**. You can call it as follows:
@@ -43,6 +46,12 @@ For example, to download the WorldPop data for Austria, you can call:
 
 In this case, the ISO 3166-1 alpha-3 must be provided for the ``<country>``.
 
+By default, the UN-adjusted 1km population dataset is used. However, you can also download the higher resolution UN-adjusted 100m population dataset as follows:
+
+    python download_worldpop.py AUT 100m
+
+### Disaggregating exposure data
+
 With WorldPop data downloaded, you can execute the core script as follows:
 
     python main_script.py <country> <group>

_config.py

@@ -24,13 +24,15 @@
 import rasterio.transform
 import rasterio.mask
 from rasterio import Affine
-from osgeo import gdal, gdalconst
 
 # Shapefiles - working with vector data
 import shapely.geometry
 from shapely.ops import cascaded_union
 import geopandas as gpd
 
+# GDAL
+from osgeo import gdal, gdalconst
+
 # Numeric
 import pandas as pd
 import numpy as np
@@ -53,8 +55,13 @@
 # data
 res = 0.05
 
-# Desired coordinate reference system
+# Desired coordinate reference system corresponding to input files
 desired_crs = "EPSG:4326"
+# Desired coordinate reference system for calculating areas (you may need to
+# look up which flat projection would work well for your region); however it's
+# not particularly important to the calculation and therefore OK to be a bit
+# inaccurate
+area_crs = "EPSG:3035" # Suitable for EU
 
 # Input shapefile parameters - field_name needs to exist in both the input
 # shapefile and the input exposure CSV file. Directly replace field_name

calcs/sampling.py

@@ -65,7 +65,7 @@ def resample_assets(df, assets, bound_names, mapped_field):
                 sampled_jdx = random.choices(i_jdx, weights=p_jdx, k=n_samples)
 
                 # Arrange into dataframe format
-                new_samples = df_jdx.iloc[sampled_jdx].copy()
+                new_samples = df_jdx.loc[sampled_jdx].copy()
                 new_samples["weights"] = w_samples
                 new_samples["taxonomy"] = i  # FIXME: should be automatic from retain_tags
 

download_worldpop.py

@@ -4,15 +4,18 @@
 
 
 # Parse function with grid; FIXME: Support other resolutions from WP beyond 1km
-def download_worldpop(iso, year=2020, res=1):
+def download_worldpop(iso, year=2020, res="1km"):
     ''' This function downloads a raster file of population estimates from
     WorldPop. At the moment, this will download the 2020 population estimate
     at a 1km resolution for the desired country. This can be later expanded to
     download the 100m dataset, or (where possible) a direct estimate of the
     number of buildings.'''
 
-    # Arrange url
-    url = f'ftp://ftp.worldpop.org.uk/GIS/Population/Global_2000_{year}/{year}/{iso.upper()}/{iso.lower()}_ppp_{year}.tif'
+    # Arrange url; defeault  - 1km, UN adj
+    url = f'ftp://ftp.worldpop.org.uk/GIS/Population/Global_2000_2020_1km_UNadj/{year}/{iso.upper()}/{iso.lower()}_ppp_2020_1km_Aggregated_UNadj.tif'
+    if res == "100m":
+        # If requested, use 100m, UN adj
+        url = f'ftp://ftp.worldpop.org.uk/GIS/Population/Global_2000_2020/{year}/{iso.upper()}/{iso.lower()}_ppp_{year}_UNadj.tif'
 
     # Submit request
     r = urllib.request.urlopen(url)
@@ -37,7 +40,9 @@ def download_worldpop(iso, year=2020, res=1):
     # Parse arguments (group name)
     iso = sys.argv[1]
     year = 2020
-    res = 1
+    res = "1km"
+    if len(sys.argv) > 2:
+        res = sys.argv[2]
 
     # Call function
     download_worldpop(iso, year, res)

main_script.py

@@ -25,22 +25,41 @@ def main(mapped_field, desired_level, country_name, country_iso, group):
     anticipated file names/directories, the desired resolution, and the
     desired coordinate reference system.'''
 
+    # --------------------------------------------------------------------------
+    #   PRINT INFO FOR USER
+    # --------------------------------------------------------------------------
+
+    # Arrange full WP path
+    wp_name = f"{country_iso.lower()}_ppp_{worldpop_year}.tif"
+    wp_path = os.path.join(wp_directory, wp_name)
+
+    # Arrange full shapefile path
+    shp_path = os.path.join(shp_directory, shp_file)
+
+    # Arrange full exposure
+    exp_name = f"Exposure_{group}_{country}.csv"
+    exp_path = os.path.join(exp_directory, exp_name)
+
+    # Print information
+    print(f"\nThis code is using the following files:")
+    print(f"- Original exposure: {exp_path}")
+    print(f"- Corresponding admin divisions: {shp_path}")
+    print(f"- WorldPop data: {wp_path}")
+    print(f"Both exposure and admin division files expected to have field {mapped_field}. Target resolution is {res}.\n")
+    print("If any of this information seems incorrect, please check input parameters within _config.py\n")
+
     # --------------------------------------------------------------------------
     #   PARSE ADMIN BOUNDS & CITIES
     # --------------------------------------------------------------------------
 
     # Read admin bounds, which must exist locally
     # This will also associate desired admin shp and associated field
-    adm = parse_adm()
+    adm = parse_adm(shp_path)
 
     # --------------------------------------------------------------------------
     #   MERGE RASTERS AND ASSOCIATE POINTS TO ADMIN BOUNDS
     # --------------------------------------------------------------------------
 
-    # Arrange full path
-    wp_name = f"{country_iso.lower()}_ppp_{worldpop_year}.tif"
-    wp_path = os.path.join(wp_directory, wp_name)
-
     # Resample (aggregate) WorldPop grid to specified coarser resolution
     wp_path = resample_raster_to_resolution(
         wp_path, wp_name, res
@@ -87,7 +106,7 @@ def main(mapped_field, desired_level, country_name, country_iso, group):
     if important_exceptions:
         print_red(f"IMPORTANT WARNING: Will not be able to properly distribute {important_exceptions}; using nearest raster values instead")
         # Add important exceptions to df and raster values from nearest point
-        df = add_excepted_bounds(df, adm[desired_level], important_exceptions,
+        df = add_excepted_bounds(df, adm, important_exceptions,
                                  mapped_field)
         # Update bound_names accordingly
         bound_names = df[mapped_field].unique()