<jupyter_start><jupyter_text>Identify settlements This notebook identifies settlements that can be used as origins in an accessibility analysis for Khbyer Pakhtunkhwa. Several different identification criteria are explored<jupyter_code>import os, sys, time

import pandas as pd
import geopandas as gpd

import shapely
from shapely.geometry import shape
from shapely.wkt import loads
from shapely.errors import TopologicalError

import rasterio
from rasterio import features, transform
from rasterio.mask import mask
from rasterio.transform import Affine
from rasterio.warp import calculate_default_transform, reproject, Resampling

import numpy as np

from rasterstats import zonal_stats

sys.path.append('../../src/')

import skimage.graph as graph
import gostrocks.src.GOSTRocks as gr
from gostrocks.src.GOSTRocks.misc import tPrint
from GOSTNets_Raster.src.GOSTNets_Raster import market_access as ma<jupyter_output><empty_output><jupyter_text>File paths<jupyter_code>geo_pth = r"/Users/robert/Documents/Jobs/WB/Pakistan/Data"
pop_pth = r"/Users/robert/Documents/Jobs/WB/Pakistan/Data/Population"
data_dir = r'../../data'
vect_dir = r'vect_inputs'
rast_dir = r'rast_inputs'
out_pth = r'../../data/vect_out'<jupyter_output><empty_output><jupyter_text>Projections<jupyter_code>dest_crs = 'EPSG:32642' # this is a Pakistani UTM projection, assign correct projection for project area<jupyter_output><empty_output><jupyter_text>Functions<jupyter_code># Lightly adapted from https://gis.stackexchange.com/questions/290030/what-does-it-mean-to-reproject-a-satellite-image-from-utm-zone-13n-to-wgs84

def reproject_tif(source_file, destination_file,dest_crs):
    """Re-projects tif at source file to destination CRS at destination file.

    Args:
        source_file: file to re-project
        destination_file: file to store re-projection

    Returns:
        destination_file: where the re-projected file is saved at
    """

    with rasterio.open(source_file) as src:
        dst_crs = dest_crs
        transform, width, height = calculate_default_transform(
            src.crs,
            dst_crs,
            src.width,
            src.height,
            *src.bounds
        )

        kwargs = src.meta.copy()
        kwargs.update({
            'crs': dst_crs,
            'transform': transform,
            'width': width,
            'height': height,
            "compress":'LZW'
        })

        with rasterio.open(destination_file, 'w', **kwargs) as dst:
            for i in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=transform,
                    dst_crs=dst_crs,
                    resampling=Resampling.nearest,
                    num_threads=-1
                )

        return destination_file<jupyter_output><empty_output><jupyter_text>Load in raster data WSF Files<jupyter_code>wsf_raw_fil = 'KP_WSF2019_binary.tif'
wsf_3x3_fil = 'KP_WSF2019_binary_skimage_modal_3by3.tif'
wsf_5x5_fil = 'KP_WSF2019_binary_skimage_modal_5by5.tif'
wsf_7x7_fil = 'KP_WSF2019_binary_skimage_modal_7by7.tif'
reproject_tif(os.path.join(data_dir,rast_dir,wsf_raw_fil),
             os.path.join(data_dir,rast_dir,'KP_WSF2019_binary_32642.tif'),
             'EPSG:32642')
with rasterio.open(os.path.join(data_dir,rast_dir,'KP_WSF2019_binary_32642.tif'),'r') as wsf_raw_raw:
    wsf_raw = wsf_raw_raw.read(1)
    wsf_raw_prof = wsf_raw_raw.profile<jupyter_output><empty_output><jupyter_text>Now WorldPop<jupyter_code>wp_kp_fil = 'Population/WorldPop/KP_Analysis_WorldPop_2020.tif'
with rasterio.open(os.path.join(geo_pth,wp_kp_fil),'r') as wp_kp_raw:
    wp_kp = wp_kp_raw.read(1)
    wk_kp_prof = wp_kp_raw.profile<jupyter_output><empty_output><jupyter_text>Settlements<jupyter_code>settle = gpd.read_file(os.path.join(data_dir,vect_dir,'KP_Settlements_NGA_2017.gpkg'),driver="GPKG")
settle = settle.set_crs(4326).to_crs(dest_crs)
# GeoNames contains duplicate values for places with alternate names (curiously). We need to drop them
settle = settle[settle.NT == 'N']
filter_list = ['PPL', 'PPLA', 'PPLA2']
settle = settle[settle['DSG'].isin(filter_list) == True]
len(settle)
settle = settle.reset_index().drop('index',axis=1)
settle.head()<jupyter_output><empty_output><jupyter_text>Transform settlements<jupyter_code>set200buff = settle.copy()
set200buff.geometry = set200buff.geometry.buffer(200)
set200buff.to_file(os.path.join(data_dir,vect_dir,"KP_Settlements_NGA_2017_200mBuff_32642.shp"))
settle.head(2)<jupyter_output><empty_output><jupyter_text>Aggregate WSF footprints within settlements<jupyter_code># run zonal stats on just the geometry column of set200, then join it back to the original GDF based on common indices

set200_zs = set200buff.join(
    pd.DataFrame(
        zonal_stats(
                 vectors = set200buff[['geometry','FID_Settle']],\
                 raster = wsf_raw,\
                 affine = wsf_raw_raw.transform,\
                 stats = 'sum'
        )
    ),
    how='left',on='FID_Settle'
)
# run zonal stats on just the geometry column of set200, then join it back to the original GDF based on common indices

settle_pts_zs = settle.join(
    pd.DataFrame(
        zonal_stats(
                 vectors = set200buff['geometry'],\
                 raster = wsf_raw,\
                 affine = wsf_raw_raw.transform,\
                 stats = 'sum'
        )
    ),
    how='left'
)
set200_zs.tail()
settle_pts_zs.tail()<jupyter_output><empty_output><jupyter_text>Filter down settlements<jupyter_code>set_any = settle_pts_zs[settle_pts_zs['sum'] > 0]
set100 = settle_pts_zs[settle_pts_zs['sum'] > 100]
set200 = settle_pts_zs[settle_pts_zs['sum'] > 200]
set500 = settle_pts_zs[settle_pts_zs['sum'] > 500]
set1000 = settle_pts_zs[settle_pts_zs['sum'] > 1000]
# 1500 and above yield no points
print(len(set_any))
print(len(set100))
print(len(set200))
print(len(set500))
print(len(set1000))<jupyter_output>10196
3863
2482
749
66<jupyter_text>Export<jupyter_code>set200_zs.to_file(os.path.join(data_dir,vect_dir,"KP_Settlements_NGA_2017_200mBuff_ZS_Results_v2.gpkg"),\
                  layer='Buffered_Polygons',driver="GPKG")

# settle_pts_zs.to_file(os.path.join(data_dir,vect_dir,"KP_Settlements_NGA_2017_200mBuff_ZS_Results.gpkg"),\
#                       layer='Points',driver="GPKG")
set_any.to_file(os.path.join(data_dir,vect_dir,"KP_NGA_Settlements_200mBuff_Filtered.gpkg"),layer="Any",driver="GPKG")
set100.to_file(os.path.join(data_dir,vect_dir,"KP_NGA_Settlements_200mBuff_Filtered.gpkg"),layer="100plus",driver="GPKG")
set200.to_file(os.path.join(data_dir,vect_dir,"KP_NGA_Settlements_200mBuff_Filtered.gpkg"),layer="200plus",driver="GPKG")
set500.to_file(os.path.join(data_dir,vect_dir,"KP_NGA_Settlements_200mBuff_Filtered.gpkg"),layer="500plus",driver="GPKG")
set1000.to_file(os.path.join(data_dir,vect_dir,"KP_NGA_Settlements_200mBuff_Filtered.gpkg"),layer="1000plus",driver="GPKG")<jupyter_output><empty_output><jupyter_text>Develop settlement catchments and aggregate populations within them We need to assign populations to settlements to capture their relative importance in the gravity model. This could be done with some sort of buffer but this is crude and would also create double counting issues when settlements are close, unless tedious and somewhat arbitrary difference calculations are imposed on the buffer polygons.A better solution is to generate exclusive catchments for the settlements based on the friction surface itself, then summarize populations of interest within these catchments. This has the double benefit of avoiding double counting and also better reflecting how market structures actually work: farmers bring produce to the local town and sell it to middlemen who transport it onwards to the main markets.We use GOST's market catchment code in combination with our friction surface to prepare the catchments Look at our settlements layer<jupyter_code>set200.head()
set100.columns
fric_dry = os.path.join(geo_pth,'Accessibility/KP_dry_friction1_210817.tif')<jupyter_output><empty_output><jupyter_text>Create an mcp object from the friction surface raster<jupyter_code>inR = rasterio.open(fric_dry)
inD = inR.read()[0,:,:] * inR.transform.a # important to specify pixel size in meters here in oder to get correct measurements
inD = np.nan_to_num(inD)
mcp = graph.MCP_Geometric(inD)<jupyter_output><empty_output><jupyter_text>Find the MCP destination indices on the raster surface and stash them in a list of tuples<jupyter_code>dests = ma.get_mcp_dests(inR, set100, makeset=False)<jupyter_output><empty_output><jupyter_text>use the dests and mcp object to generate catchment areas around the settlements<jupyter_code>ma.generate_market_sheds??<jupyter_output><empty_output><jupyter_text>Generate the actual catchments -- this will take a while<jupyter_code>set100.index
set100_catch = ma.generate_market_sheds(inR,set100,inR.transform.a)
with rasterio.open(geo_pth,'Accessibility/KP_Set200_Marketsheds_test.tif',**inR.profile) as dst:
    dst.write(set100_catch)
# set100_catch = ma.generate_feature_vectors(inR,mcp,set100,[30],'FID_Settle') # veeery slow<jupyter_output>14:10:43	1 of 3863: 44285
14:30:14	2 of 3863: 44199<jupyter_text>OLD Join the catchments to the original dataset<jupyter_code>set100_catch.head()
set100_catch_fin = gpd.GeoDataFrame(pd.merge(set100_catch,set100.drop('geometry',axis=1),how='left',on='FID_Settle'))
# run zonal stats on just the geometry column of set100, then join it back to the original GDF based on common indices

set100_catch_zs = set100.join(
    pd.DataFrame(
        zonal_stats(
                 vectors = set100_catch,\
                 raster = wp_kp,\
                 affine = wp_kp_raw.transform,\
                 stats = 'sum'
        )
    ),
    how='left',on='FID_Settle'
)<jupyter_output>/Users/robert/.pyenv/versions/geo/lib/python3.9/site-packages/rasterstats/io.py:302: UserWarning: Setting nodata to -999; specify nodata explicitly
  warnings.warn("Setting nodata to -999; specify nodata explicitly")
{'fraction_non_alphanumeric': 0.09945022580011781, 'fraction_numerical': 0.035637149028077755, 'mean_word_length': 4.521409214092141, 'pattern_counts': {'":': 1, '<': 66, '<?xml version=': 0, '>': 71, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '363961', 'license': 'CC0-1.0', 'n_tokens_mistral': 3828, 'n_tokens_neox': 3480, 'n_words': 692, 'path': 'notebooks/Data_Prep/Data Prep - Identify Settlements.ipynb', 'repo_name': 'worldbank/Khyber-Pakhtunkhwa-Accessibility-Analysis', 'size': '10186'}