Notebook

In [10]:

import os
import glob
import csv
import requests
from pathlib import Path
import xarray as xr
import rioxarray as rio
import pandas as pd
import geopandas as gpd
import numpy as np
import rasterio
import folium
from rasterio.features import shapes
from branca.colormap import linear

In [2]:

city = 'Luanda'

In [3]:

data_folder = Path(r"D:\World Bank\CRP\data\WaterGAP")
output_folder = Path('output')
html_folder = Path('html')
int_output_folder = Path('intermediate_output')

In [4]:

# read catchment AOI
aoi = gpd.read_file('AOI/luanda_catchment_level4.shp').to_crs(epsg = 4326)

In [14]:

var_list = {'ncrun': 'Net cell runoff', 
            'qs': 'Fast surface and fast subsurface runoff'}
stat_list = {'mean': 'mean',
             'std': 'standard deviation'}
year_list = range(1987, 2017)

Convert netcdf to raster¶

In [12]:

for var in var_list:
    if not os.path.exists(data_folder / f'{var}_raster2d'):
        os.mkdir(data_folder / f'{var}_raster2d')
    
    # open netcdf
    nc = rio.open_rasterio(data_folder / f'watergap_22d_WFDEI-GPCC_histsoc_{var}_monthly_1901_2016.nc4',
                           decode_times = False)
    # set time dimension
    units, reference_date = nc.time.attrs['units'].split('since')
    nc['time'] = pd.date_range(start = reference_date, periods = nc.sizes['time'], freq = 'MS')
    
    # output raster 2D if it does not already exist
    for year in year_list:
        for month in range(1, 13):
            output_raster = data_folder / f'{var}_raster2d' / f'{var}_{str(year)}{str(month).zfill(2)}.tif'
            if not os.path.exists(output_raster):
                nc.rio.write_crs("epsg:4326", inplace = True).sel(time = f'{str(year)}-{str(month).zfill(2)}-01').rio.to_raster(output_raster)

In [40]:

# clip each raster
for var in var_list:
    for year in year_list:
        for month in range(1, 13):
            input_raster = data_folder / f'{var}_raster2d' / f'{var}_{str(year)}{str(month).zfill(2)}.tif'
            with rasterio.open(input_raster) as src:
                # shapely presumes all operations on two or more features exist in the same Cartesian plane.
                out_image, out_transform = rasterio.mask.mask(
                    src, aoi.geometry, crop = True, all_touched = True)
                if year == year_list[0] and month == 1:
                    out_meta = src.meta.copy()

            out_meta.update({"driver": "GTiff",
                             "height": out_image.shape[1],
                             "width": out_image.shape[2],
                             "transform": out_transform})

            output_file = f'{var}_{str(year)}{str(month).zfill(2)}.tif'
            with rasterio.open(int_output_folder / output_file, "w", **out_meta) as dest:
                dest.write(out_image)

In [47]:

# average raster (or other summary stats) and save
for var in var_list:
    raster_list = []
    
    for year in year_list:
        for month in range(1, 13):
            with rasterio.open(int_output_folder / f'{var}_{str(year)}{str(month).zfill(2)}.tif') as src:
                array = src.read(1)
                array[array == src.meta['nodata']] = np.nan
                raster_list.append(array)
                if year == year_list[0] and month == 1:
                    meta = src.meta
                    meta['nodata'] = np.nan
    
    output_raster_mean = np.nanmean(raster_list, axis = 0)
    output_raster_std = np.nanstd(raster_list, axis = 0)
    
    with rasterio.open(output_folder / f'{city.lower()}_{var}_mean.tif', 'w', **meta) as dst:
        dst.write(output_raster_mean, 1)
    with rasterio.open(output_folder / f'{city.lower()}_{var}_std.tif', 'w', **meta) as dst:
        dst.write(output_raster_std, 1)

C:\Users\Owner\AppData\Local\Temp\ipykernel_11236\2328891307.py:15: RuntimeWarning: Mean of empty slice
  output_raster_mean = np.nanmean(raster_list, axis = 0)
C:\Users\Owner\anaconda3\envs\crp\lib\site-packages\numpy\lib\nanfunctions.py:1879: RuntimeWarning: Degrees of freedom <= 0 for slice.
  var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof,

In [42]:

# generate time series stats
for var in var_list:
    time_series = {'mean': {}, 'std': {},
                   'max' : {}, 'min': {}}
    
    for year in year_list:
        for month in range(1, 13):
            with rasterio.open(int_output_folder / f'{var}_{str(year)}{str(month).zfill(2)}.tif') as src:
                src_array = src.read(1)
                src_array = src_array[src_array != src.meta['nodata']]
                time_series['mean'][f'{str(year)}{str(month).zfill(2)}'] = np.nanmean(src_array)
                time_series['std'][f'{str(year)}{str(month).zfill(2)}'] = np.nanstd(src_array)
                time_series['min'][f'{str(year)}{str(month).zfill(2)}'] = np.nanmin(src_array)
                time_series['max'][f'{str(year)}{str(month).zfill(2)}'] = np.nanmax(src_array)

    with open(output_folder / f'{city.lower()}_{var}.csv', 'w') as f:
        f.write('year,month,mean,std,min,max\n')
        for year in year_list:
            for month in range(1, 13):
                f.write('%s,%s,%s,%s,%s,%s\n' % (str(year), str(month).zfill(2),
                                                 time_series['mean'][f'{str(year)}{str(month).zfill(2)}'],
                                                 time_series['std'][f'{str(year)}{str(month).zfill(2)}'],
                                                 time_series['min'][f'{str(year)}{str(month).zfill(2)}'],
                                                 time_series['max'][f'{str(year)}{str(month).zfill(2)}']))

In [10]:

# delete intermediate outputs
for var in var_list:
    for year in year_list:
        for month in range(1, 13):
            os.remove(int_output_folder / f'{var}_{str(year)}{str(month).zfill(2)}.tif')

Data conversion¶

In [48]:

for var in var_list:
    for stat in stat_list:
        with rasterio.open(output_folder / f'{city.lower()}_{var}_{stat}.tif') as src:
            image = np.float32(src.read(1))
            results = ({'properties': {'raster_val': np.float32(v)}, 'geometry': s} for i, (s, v) in enumerate(shapes(image, mask = None, transform = src.transform)))
            geoms = list(results)
            gpd_polygonized_raster = gpd.GeoDataFrame.from_features(geoms)
            gpd_polygonized_raster = gpd_polygonized_raster[gpd_polygonized_raster.raster_val != 0]

        gpd_polygonized_raster.to_file(output_folder / f'{city.lower()}_{var}_{stat}.geojson', driver = 'GeoJSON', crs = 'EPSG:4326')

Map¶

In [88]:

# Create a folium map centered on the AOI
mymap = folium.Map(control_scale = True)
aoi_bounds = aoi.total_bounds.tolist()
mymap.fit_bounds([aoi_bounds[:2][::-1], aoi_bounds[2:][::-1]])

folium.GeoJson(
    aoi,
    style_function=lambda feature: {
        'fillColor': 'transparent',
        'color': 'gray',
        'weight': 3
    },
    control = False).add_to(mymap)

Out[88]:

<folium.features.GeoJson at 0x1e8eac5d3c0>

In [89]:

for var in var_list:
    for stat in stat_list:
        # read json
        json = gpd.read_file(output_folder / f'{city.lower()}_{var}_{stat}.geojson').to_crs('epsg:4326')
        
        # set colormap
        if var == 'ncrun':
            colormap = linear.GnBu_09.scale(json['raster_val'].min(), json['raster_val'].max())
        elif var == 'qs':
            colormap = linear.YlGnBu_09.scale(json['raster_val'].min(), json['raster_val'].max())
        
        # determine whether to show layer by default
        show = False
        if var == 'ncrun' and stat == 'mean':
            show = True
        
        # Add GeoJson layer with choropleth style
        folium.GeoJson(
            json,
            name = f'{var_list[var]} ({stat_list[stat]})',
            # overlay = False,
            style_function=lambda feature: {
                'fillColor': 'transparent' if pd.isna(feature['properties']['raster_val']) else colormap(feature['properties']['raster_val']),
                # 'fillColor': colormap(feature['properties']['raster_val']),
                'fillOpacity': 0.7,
                'weight': 0,
            },
            highlight_function=lambda x: {'fillOpacity': 0.9},
            tooltip=folium.GeoJsonTooltip(fields=['raster_val'], aliases = [f'{var_list[var]} ({stat_list[stat]})'], 
                                          labels=True, sticky=False),
            show = show
        ).add_to(mymap)

# Add Layer Control to toggle the choropleth layer
folium.LayerControl().add_to(mymap)

Out[89]:

<folium.map.LayerControl at 0x1e8ea8fa860>

In [90]:

mymap

Out[90]:

Make this Notebook Trusted to load map: File -> Trust Notebook

In [91]:

mymap.save(html_folder / "runoff.html")