Olivia Burge review

docs/methodology.md

@@ -49,7 +49,9 @@ We used the land cover information from the Landcover Database version 5.0 for t
 The data is downloaded as a vector, rasterised as a binary raster at 10 m² resolution, and then downsampled to 100 m² with a summation algorithm to obtain a 0–100 value for measuring partial pixel cover at this scale (a value of 20 indicates that 20% of the pixel is covered in cropland or forestry). These values are converted to a footprint score using the equation:
 
 $$
-F = \frac{x}{100}\cdot4
+F = \begin{cases} 7 & \text{if $x > 20$} \\
+4 & \text{if $0 < x \le 20$} \\
+0 & \text{otherwise} \end{cases}
 $$
 
 ### Pasture
@@ -59,9 +61,7 @@ Pasture data is obtained similarly to croplands. Classes 2 (urban parkland/open
 The cover values are obtained in the same fashion as cropland and are converted to a footprint score using the rule:
 
 $$
-F = \begin{cases} 7 & \text{if $x > 20$} \\
-4 & \text{if $0 < x \le 20$} \\
-0 & \text{otherwise} \end{cases}
+F = \frac{x}{100}\cdot4
 $$
 
 ### Roads
@@ -131,7 +131,11 @@ $$
 
 ### Final human footprint index
 
-All eight components are combined additvely. The New Zealand coastline (1:50,000 scale) is re-used in this final step in order to firmly set no-data values for marine spaces. Outside of marine areas, the minimum value is 0. This results in a raster layer where all terrestrial spaces have a minimum value of 0 (wilderness) up to a hypothetical maximum of 61. The datatype is 32-bit floating-point, rather than an integer, because some components use logarithmic or exponential functions and the values are not rounded.
+All eight components are combined additvely.
+
+The New Zealand coastline (1:50,000 scale) is re-used in this final step in order to firmly set no-data values for marine spaces. Outside of marine areas, the minimum value is 0. This results in a raster layer where all terrestrial spaces have a minimum value of 0 (wilderness) up to a hypothetical maximum of 50. Following [[19]](#19), urban, cropland, and pasture follow a "land-use exclusion principle" whereby the three types may not co-occur. Although the individual footprint values are calculated independently, at the combination stage the first non-zero value among urban, cropland, and wilderness footprint components is taken, in that order. (If this were not the case, and summation was performed naïvely, the hypothetical maximum value would be 61.)
+
+The output datatype is 32-bit floating-point, rather than an integer, because some components use logarithmic or exponetnial functions and the values are not rounded.
 
 ## Reproducibility
 

images/HFI-2018.pgw

@@ -1,6 +1,6 @@
-591.7775015612161269
+606.55732868093252819
 0
 0
--591.7775015612161269
-2238589.48885078076273203
-7748253.92202259693294764
+-606.55732868093252819
+2017329.46796434046700597
+7755937.42878835927695036

src/Snakefile

@@ -22,7 +22,7 @@ NZ_COAST_RASTER = OUTD / 'data/downloads/nz-coastlines-and-islands-polygons-topo
 TARGET_YEARS = ['2012', '2018']
 
 FOOTPRINT = OUTD / 'data/footprints/output/human_footprint_nz-{year}.tif'
-ARCHIVE = OUTD / 'archival/human_footrint.archive.{year}.zip'
+ARCHIVE = OUTD / 'archival/human_footrint.archive.zip'
 
 include: "rules/lcdb.smk"
 include: "rules/worldpop.smk"
@@ -42,14 +42,15 @@ rule all:
 source=lambda wc, input: pathlib.PurePosixPath(input.source).relative_to(OUTD),
 
 rule archive:
-    input: FOOTPRINT, GLOBAL_HFP_DIFF
+    input: expand(FOOTPRINT, year=TARGET_YEARS), expand(GLOBAL_HFP_DIFF, year=TARGET_YEARS)
     output: ARCHIVE
     params:
         target=OUTD / 'data'
     shell: '''
         mkdir -p $(dirname {output})
         zip -r {output} {params.target}
     '''
+
 rule download_nz_coastlines_and_islands_polygon:
     output: NZ_COAST_VECTOR
     threads: 2
@@ -89,15 +90,16 @@ rule rasterise_nz_coastlines_and_islands_polygon:
         '''
 
 rule summation:
-    input: NZ_COAST_RASTER, WORLDPOP_FOOTPRINT_NZ, GHS_FOOTPRINT, VNL_FOOTPRINT, PASTURE_FOOTPRINT, CROPLAND_FOOTPRINT, ROADS_FOOTPRINT, RAIL_FOOTPRINT, NAVIGABLE_WATER_FOOTPRINT
+    input: NZ_COAST_RASTER, WORLDPOP_FOOTPRINT_NZ, VNL_FOOTPRINT, GHS_FOOTPRINT, CROPLAND_FOOTPRINT, PASTURE_FOOTPRINT, ROADS_FOOTPRINT, RAIL_FOOTPRINT, NAVIGABLE_WATER_FOOTPRINT
     output: FOOTPRINT
     conda: './envs/gdal.yml'
     log: LOGD / 'summation-{year}.log'
     params:
         tmp=TMPD,
-        max_value=10+10+10+7+4+8+8+4, # Maximum hypothetical value
+        max_value=50, # Maximum hypothetical value
         input_layers_calc=lambda wildcards, input: ' '.join(map(lambda x: f'-{x[0]} {x[1]}.vrt', zip(list(string.ascii_uppercase[:len(input)]), input))),
-        input_calc=lambda wildcards, input: '+'.join(list(string.ascii_uppercase)[1:len(input)]),
+        # input_calc=lambda wildcards, input: '+'.join(list(string.ascii_uppercase)[1:len(input)]),
+        input_calc='B+C+((D*(D>0))|(E*(E>0))|(F*(F>0)))+G+H+I',
         creation_options=" ".join(f'--co {k}={v}' for k, v in config['compression_co']['zstd_pred3'].items())
     shell: '''
         mkdir -p $(dirname {output})
@@ -107,7 +109,7 @@ rule summation:
         gdal_calc.py {params.input_layers_calc} --calc="(A==0)*-1+(A==1)*({params.input_calc})" --extent=union \
             --outfile {params.tmp}/$(basename {output}) --overwrite \
             {params.creation_options} --NoDataValue=-1
-        gdal_calc.py -A {params.tmp}/$(basename {output}) -B {input[0]} --calc="where(logical_and(A>{params.max_value},B==1),0,A)" \
+        gdal_calc.py -A {params.tmp}/$(basename {output}) -B {input[0]} --calc="where(logical_and(A>{params.max_value},B==1),{params.max_value},A)" \
             --outfile {output} --overwrite \
             {params.creation_options} --hideNoData --NoDataValue=-1
         gdal_edit.py -stats {output}

src/rules/eog.smk

@@ -17,20 +17,23 @@ VNL_URLS = {
     2012: 'https://eogdata.mines.edu/nighttime_light/annual/v21/2012/VNL_v21_npp_201204-201212_global_vcmcfg_c202205302300.median_masked.dat.tif.gz',
 }
 
+DECILE_BASELINE_YEAR = 2012
+
 VNL_YEARS = list(map(str, VNL_URLS.keys()))
 
 # Median monthly radiance, nW/cm^2/sr
-# NB output is scaled 0-255 (Byte) over New Zealand to make computation of deciles computationally easier for the Human Footprint index, 
-# and therefore is not actually in units of nW/cm^2/sr, but the intermediate data is retained if needed.
-# Since Human Footpring mapping is using VNL data for the computation of deciles, this conversion does not result in lost information. 
+# NB output is scaled to 0-65535 (UInt16) over New Zealand to make computation of deciles computationally easier for the Human Footprint index,
+# and therefore is not actually in units of nW/cm^2/sr, but the original (clipped, floating point) data is retained if needed.
+# Since the Human Footprint index converts VNL data into deciles, this conversion does not result in lost information,
 rule download_project_clip_vnl:
     output: VNL
     wildcard_constraints:
         year='\d{4}'
     params:
         url=lambda wildcards: VNL_URLS[get_nearest(VNL_URLS, wildcards.year)],
         extent=config['extent'],
-        creation_options=" ".join(f'-co {k}={v}' for k, v in config['compression_co']['zstd_pred3'].items())
+        creation_options=" ".join(f'-co {k}={v}' for k, v in config['compression_co']['zstd_pred3'].items()),
+        scale_max=1000 # nW/cm^2/sr
     conda: '../envs/gdal.yml'
     log: LOGD / "download_project_clip_vnl_{year}.log"
     shell: '''
@@ -42,19 +45,18 @@ rule download_project_clip_vnl:
             -multi -wo NUM_THREADS=ALL_CPUS \
             {output}.4326.tif {output}.unscaled.tif
         gdal_edit.py -stats {output}.unscaled.tif
-        gdal_translate -ot Byte -scale {output}.unscaled.tif {output}
+        gdal_translate -ot Uint16 -scale 0 {params.scale_max} 0 65535 {output}.unscaled.tif {output}
         gdal_edit.py -stats {output}
     '''
 
-# NB perform "gdalinfo -hist {output}" to verify that the output is in 10 approximately equal bins (excluding 0).
-# The result won't have exactly bins: values stradling the edge aren't distributed evenly between boundary values,
-# Yet, for 2020, each class 1-10 has between 131,825 to 131,859 pixels, and the remaining 1,004,044,827 pixels are 0.
 rule footprint_vnl:
-    input: VNL
+    input:
+        night_light=VNL,
+        baseline=expand(VNL, year=[DECILE_BASELINE_YEAR])
     output: VNL_FOOTPRINT
     log: LOGD / "footprint_vnl_{year}.log"
     threads: 5
     params:
-        logLevel='DEBUG'
+        logLevel='INFO'
     conda: '../envs/rasterio.yml'
     script: '../scripts/decile.py'

src/rules/lcdb.smk

@@ -101,6 +101,6 @@ rule footprint_pasture:
         creation_options=" ".join(f'--co {k}={v}' for k, v in config['compression_co']['zstd_pred2'].items())
     shell: '''
         mkdir -p $(dirname {output})
-        gdal_calc.py --outfile={output} -A {input} --calc="A/100.0*4" --overwrite {params.creation_options}
+        gdal_calc.py --type Float32 --outfile={output} -A {input} --calc="A/100.0*4" --overwrite {params.creation_options}
         gdal_edit.py -stats -a_srs EPSG:3851 {output}
         '''

src/scripts/decile.py

@@ -1,6 +1,8 @@
 """
 Convert a numerical raster into a 1-10 scale, on the basis of ten equal quantiles.
 0 and NaN are excluded from the quantile calculation, and retained in the output.
+Quantiles can be calculated on the basis of a different raster than the one used
+for classification, e.g. a baseline year.
 """
 
 from concurrent.futures import ThreadPoolExecutor
@@ -47,20 +49,23 @@ def trim_zeros(filt: np.ndarray, trim='fb') -> np.ndarray:
     return filt[i:j]
 
 def calculate_deciles(data: np.ndarray) -> np.ndarray:
-    return np.percentile(
+    deciles = np.quantile(
         trim_zeros(np.sort(data.flatten(), axis=0, kind='stable'), trim='f'),
-        np.arange(10, 100, 10)
+        np.arange(10, 100, 10),
+        method='closest_observation'
     )
+    return np.insert(deciles, 0, 0) # Add the bottom value
 
 def main(workers=smk.threads):
 
-    with rio.open(smk.input[0]) as ds:
+    baseline = Path(smk.input['baseline'][0])
+    target = Path(smk.input['night_light'])
+    with rio.open(baseline) as ds:
         deciles : np.narray = calculate_deciles(ds.read(1, masked=False))
-    deciles = np.insert(deciles, 0, 0)
     logging.info(deciles)
 
-    logging.info(f"Opening {smk.input[0]} for reading")
-    with rio.open(smk.input[0]) as src:
+    logging.info(f"Opening {target} for reading")
+    with rio.open(target) as src:
         # Create a destination dataset based on source params. The
         # destination will be tiled, and we'll process the tiles
         # concurrently.